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ASSTMCT 


A three-dimensional, nonlinear nozzle admittance relation is developed 
"by solving the wave equation describing finite -amplitude oscillatory flow 
inside the subsonic portion of a choked, slowiy-convergent axisymmetric nozzle. 
This nonlinear nozzle admittance relation is then used as a "boundary condition 
in the analysis of nonlinear combust i*^!^ instability in a cylindrical liquid 
rocket combustor . In both nozzle and chamber analyses solutions are obtained 
using the Galerkin method with a series expansion consisting of the first 
tangential, second tangential, and first radial modes. Using Crocco s time-lag 
model to describe the distributed unsteady combustion process, combustion 
instability calculations are presented for different values of the following 
parameters: (l) time-lag, (2) interaction index, (s) steady-state Mach number 

at the nozzle entrance, and (4) chamber length-to -diameter ratio. In each case, 
limit-cycle pressure amplitudes and waveforms are shown for both linear and 
nonlinear nozzle admittance conditions. These results show that when the 
amplitudes of the second tangential and first radial modes are considerably 
smaller than the amplitude of the first tangential mode the inclusion of 
nozzle nonlinearities has no significant effect on the limiting aniplitude and 
pressure waveforms. 
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SUMmRY 


Recently, a three-dimensional, nonlinear nozzle admittance relation has 
teen developed. In this analysis, the wave equation for an axisymmetric, choked 
nozzle was solved using the Galerkin method, with an approximating series solu- 
tion for the velocity potential perturbation which was compatible with recent 
nonlinear combustion instability theories. Assuming, that the amplitude of the 
fundamental mode is considerably larger than- the amplitudes of the remaining 
modes in the series expansion, nonlinear admittance coefficients were determined 
as a function of the frequency and amplitude of the fundamental mode. 

The nonlinear nozzle theory was then applied in the analysis of nonlinear 
combustion instability in a cylindrical coaliustor with uniform injection of pro- 
pellants at one end and a slowly converging nozzle at the other end. The dis- 
tributed unsteady combustion process was described by means of Crocco's tame- 
lag model. The Galerkin method was used to determine the behavior of the pres- 
sure perturbation in the rocket conibTastor, where the nonlinear nozzle admxttance 
relation was used as the boundary condition at the nozzle end of the chamber. 

Tn these computations, a three-mode series expansion consisting of the first 
tangential (IT), second tangential' (2T) , and first radial (1R) modes was used. 
Since the amplitude and frequency of the Id? mode upon which the nonlinear nozzle 
admittances depend are not known a priori, an iterative solution technique 
was used. 

Combustion instability calculations have been made for different values 
of the following parameters: (l) time-lag, (2) interaction index, ( 3 ) steady 
state Mach number at the nozzle entrance, and (4) chamber length-to-diameter 
ratio. In each case limit-cycle pressure amplitudes and waveforms were 
obtained with both the linear and nonlinear nozzle admittances. These results 
show that under the' assumptions of the analysis the effect of nozzle non- 
linearities can be safely neglected in nonlinear stability calculations. 



IDJTRODUCriOM 


Various aerospace propulsion de-vices 3 such as liquid and solid propellant 
rocket motors and air hreathing jet engines, are often subject to combustion 
instabilities which are detrimental to the performance and safety of operation 
of these devices. In order to design stable engines, capabilities for a 
priori determination of the linear and nonlinear characteristics of the 
instability and the range of operating conditions for -which these engines 
are dynamically stable must be acquired. In order to perform such an analysis, 
the behavior of the e2cb.aust nozzle under oscillatory flow conditions must be 
xmderstood. In particular, it is necessary to know how a -wave generated in 
the combustion chamber is partially transmitted and partially reflected at the 
nozzle entrance. The information is usually eispressed as a boundary condition 
(usually referred to as a Ifozale Admittance Belation) that must be satisfied 
at the nozzle entrance. 

Before such a boundary condition can be derived, the nature of the -wave 
motion inside the nozzle must be investigated. The behavior of oscillations 
in a converging-diverging supercritical nozzle -was first treated by Tsien^ 
who considered the case in which the oscillation of the incoming flow is one- 
dimensional and isothermal. Crocco extended Tsien’s work to cover the 
more general cases of non- isothermal one- and three-dimensional oscillations. 
The analyses of Tsien and Crocco are both restricted to small-amplitude 
(i.e., linear) oscillations. More recently, a nonlinear nozzle theory has 
been developed by Zinn and Crocco who extended the previous linear 

theories to the investigation of the behavior of finite -amplitude waves. 

In recent st-udies conducted by Zinn, Powell, and Lores, theories were 

developed which describe the nonlinear behavior of longitudinal'^^® and 

trans-verse instabilities in liquid -propellant rocket chambers with quasi- 

steady nozzles . These theories have now been extended to situations In which 

the instabilities are three-dimensional and the rocket combustors -are attached 

11 

to conventional nozzles . All of these theories have successfully predicted 
the transient behavior, nonlinear waveforms, and limit-cycle amplitudes of 
longitudinal and tangential instabilities in unstable motors. 

In order to assess the miportance of nozzle nonlinear it ies upon the 
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aanlinear stability characteristics of various propulsion devices, a new non- 
linear nozzle theory is needed for the following reasons, Pirst, the nonlinear 
analysis of Zinn"^’ is mathematically complicated and requires considerable 
conputer time. For this reason, Zinn^s analysis has never been used to per- 
fom actual computations of the wave structure in the nozzle or the nonlinear 
nozzle response. Secondly, the nonlinear nozzle admittance .relation developed 
by Zinn is not compatible with the recently developed nonlinear combustion 
theories (see References 7 through ll) , Consequently, a linear nozzle 
bomdary condition or a short nozzle (quasi-steady) assumption had to be used 
in all of the nonlinear combustion instability theories developed to date. 

The use of a linear nozzle boundary condition in these nonlinear theories was 
justified by assuming that under the conditions of moderate amplitude oscilla- 
tions and small mean flow Mach number the effect of nozzle nonlinearities is of 
higher order and can be neglected. Thus a nonlinear nozzle analysis is needed 
to determine the validity of this assumption. Furthermore, in the case of 
transverse instabilities the "linear" nozzle has been toiown to exert a 
destabilizing effect; in these cases it is especially iuoportant to know how 
nonlinearities affect the nozzle behavior. 

Thus a nonlinear nozzle admittance relation has been developed and has 
been applied as a boundary condition in the recently-developed nonlinear 
combustion instability theories. The development of this theory, its 
application in the chamber stability analysis, and typical results for liquid- 
propellant rockets will be described in the following sections. 
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SIMBOIS 

axially dependent amplitude functions in Eq. {h) 

time dependent amplitude functions in Eq. (18) 

nozzle boundary residiial (see Eq. (lO)) 
complex ao ci al acoustic eigenvalue 

dimensionless sonic velocity, c*/c* 
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E (!') 


m 


m 

n 

P 


residual of Eq. (2) 
residual of Eq. ( 17 ) 
imaginary unit, V-1 

Bessel function of the first kind, order m 

multiple of fundamental frequency 

azimuthal mode number 
pressxire interaction index 
ataansloBless pressure, YP*/p!=f 


o o 


* , * 


f 

dimensionless radial coordinate, r /r 


* 


chamber radius 


mn 


u 


dimensionless transverse mode acoustic frequency 

t 


dimensionless time, 


d o 


* / ■J5' 


dimensionless axial velocity, u /c 


o 


ijli 

linear admittance for the p mode 


* , * 


dimensionless axial coordinate, z /r 


specific heat ratio 


th 


nonlinear admittance for the p mode 


linear admittance function 
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P 


azimuthal coordinate 

... 

damensaonless density, p /p^ 


* 


(//c*) 

d o 


T 


dimensionless pressure sensitive time lag, 
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i 

If 

U) 


steady state potential function 
Telocity, potential 
steady state stream, function 
dimensionless frequency 


Subscripts : 

e 

n 

r, i 
•ff 

o 

jZ,t 


eYaluated at the nozzle entrance 
radial mode nuaiber 

real and. imaginary parts of a con^lex quantity, respective 3;; 
evaluated at the nozzle "wall 
stagnation quantity 

partial differentiation with respect to cp,ijr,r,0 ,z, or t, 
respectively 


Superscripts : 

( )' 

o 


perturbation quantity 
steady state quantity 
dimensional quantity, complex conjugate 



approximate solution 


UOZZLE MAEfSIS 


The development of the nonlinear nozzle theory is described in detail 
in Hefs. (12) and, (13), therefore only a" brief summary will he- given in 
this section. 


5 



Development of* the Nozzle Wave EqTiation. 

As in the Zinn-Crocco analysis, ^ finite -amplitude, periodic oscilla- 
tions were assumed to occur inside the slowly convergent, subsonic portion of 
an axisymmetric nozzle operating in the supercritical range. Ihe flow in the 
nozzle was assumed to he adiahatic and inviscid and to have no body forces 
or chemical reactions. The fluid was also assumed to be calorically perfect. 
Under the further assumption of isentropic and irrotational flow the continuity 
and momentum equations were combined to obtain the following equation which 
describes the behavior of the velocity potential: 

- ~ + (y-1) 

+ (V§’V$) V^§ + i v§*v(vi*v$) 


These equations are consistent with those used in the second-order nonlinear 
combustion instability theory developed by Powell, Zinn, and Lores (see 
References 7 and lO) , 

A nozzle wave equation was obtained from Eq. (l) by expressing the 


velocity potential as the sum of a steady state and a pertxirbation 
(i.e. 5 = I + ? ), introducing the coordinate system used by Zinn 

and Crocco (see Figure l) , assuming a slowly convergent nozzle and one- 
dimensional mean flow, and neglecting third order nonlinear terms . This wave 
equation is given by: 




- 2 +f, (cp)§ 

cpfc -2 t- 

u 


tt 




■ 9 e ®0t 
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where 


+ (y+i)S^ If *09 

+ fjtf) (5')® + fg(9) t(ip^ + fg(<p) (*e')‘" 

+ (Y-i) «Pp «; - f4(<p) J ' il 
* (Y-i)|[2(tip + *;)+i*'ee]*; 

+ (Y-I)pu [2 (<.ip + «'^) + »ee] *'tp } ' ° 


f2<'P) 

^3(9) 

t> 


-2 -2 

c - -u 


-2 dcp 


u 


_ (V-<1) du^ 
2c^ dtp 


- 2 -, -2 





( 3 ) 
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figure 1. Coordinate System used for the Solution of the Oscillatory Nozzle Flow, 





Method of Solution 


In the nonlinear com'b'UBtxon instability theories deTeloped. hy Powell and 
Zinn (see Refs, 7 " H) the governing equations were solved "by means of an 
approximate solution technique known as the Galerkin Method, which is a 
special case of the Method of Weighted Eesiduals^^’^^. In these investiga- 
tions it was shown that the Galerkin Method could he successfully applied in 
the solution of nonlinear combustion instability probleias; its application 
was straightforward and it required relatively little computation time. 

Thus the Galerkin Method was also used in the nozzle analysis to determine 
the nonlinear nozzle admittance relation. 

The first step in using the Galerkin Method in the solution of the wave 
equation (i.e., Eq. (2)) was to express the velocity potential, i^, as an 
approximating series expansion. The structure of this series expansion was 
guided by the experience gained in the nonlinear nozzle admittance studies 
performed by Zinn and Crocco (see Ref, 5 ) well as in the nonlinear com- 

j 

bustion instability analyses of Powell and Zinn (see Ref. lO) . Thus the 
velocity potential was expressed as follows: 



where the , functions A (cp) are unknown coirrplex functions of the axial 
variable and 0- and ijr -dependent eigenTunctions were detemined from 

the first-order (i.e., linear) solutions by Zinn"^, For each value of the 
index p, there corresponds the mode numbers m(p) and n(p) as well as 
the number k^. This correspondence is illustrated in the table below for 
a three-term expansion consisting of the first tangential (1T) , second 
tangential (2T), and first radial (1R) modes. 
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Table 1. Three ~Mode Esjjansion 
p m(p) n(p) kp Mode 


111 


1 ^ IT 


2 2 1 


2 ■ 2T 


3 0 


1 



in 


In the t hue-dependence, u> is the fundamental frequency which must he 
specified and the integer gives the frequency of the higher harmonics. 

The values of k for the various modes appearing in Eq. ( 4 ) were determined 

ir 

from the results of the nonlinear comhustion instability analysis of Powell 
and Zinn^^. for example it was fotuid that, due to nonlinear coupling between 
modes, the 2T and IR modes oscillated with twice the frequency of the IT mode. 
Thus in Eq. ( 4 ) = 1 and k^ = ~ amplitudes and phases of the 

various inodes depend on the axial location (i.e., <p) in the nozzle through the 
unknown functions A (cp) . 

Next the assumed series expansion for (i.e., Eq. ( 4 )) was substituted 
into the wave equation (i.e., Eq. (2)) to form the residual, . According 

to the Galerkin method, the residual was required to satisfy the 

following orthogonality conditions : 



( 5 ) 


j = 1, 2, ... N... 


where N is the number of terms in the series expansions of the dependent 
variables.' The weighting ipmctions in Eq. ( 5 ) correspond to the assumed time 
and space dependences of the terms that appear in the series esqiansion. 
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The time integration is performed over one period of oscillation, T = 2tt/uj, 

•while the spatial integration is performed over any surface of 9 = constant 

in the nozzle (in Eq. (5) dS indicates an incremental area on this snrfaee) . 

Evaluating the spatial and temporal integrals in Eg,. (5) yielded a 

system of K nonlinear, second order, coupled, complex ordinary differential 

equations to be solved for the complex amplitude functions A (cp), 

E 

Unfortunately these equations were not quasi-linear; that is, the highest 

order derivatives appeared in the nonlinear terms. This greatly complicated the 

numerical solution of these equations, thus an additional approximation 'was 

made to obtain a quasi-linear system, of equations. 

This additio3oal approximation was based on the well-known fact that most 

transverse instabilities behave like the first tangential (IT) mode. Based 

11 

on the results of the recent nonlinear coaibustion instability theory , it was 
assumed that the amplitude of the IT mode was considerably larger. than the 
amplitudes of the remaining modes in the series solution. Through an order of 
magnitude analysis correct to second order, the original non-quasilinear system 
of equations was reduced to the following linear inhomogeneous system of 
equations : 


dA. 

H^(9) + M^(cp) — + H^(9)A^(cp) = 0 (6) 

dcp dcp 


d^A dA 

H (9) + M (9) — ^ + U (9)A (9) 

-y d.9 d9 




} 


p = 2, 3, ... U 


where 



(7) 


Hp(9) = u^(c^ - u^) 

-2 

M. (<P) = - i? — + 2ik col 

» acp - '^’-^ ■ 

!.,(.) = [- i . p ^ . ^y] 

and Ip are inhomogeneous tenns which are Actions of cp and the amplitude 
of the II modes * 

It can he seen that the above equations are decoupled with respect 
to the IT mode; that is, the solution for can be obtained independently 
of the amplitudes of the other modes. Thus to second order the nozzle non- 
linearities do not affect the IT mode. On the other hand, the nozzle non- 
linearities influence the aiig)litudes of the higher modes (i.e., A^ a.nd A^) 
by means of the inhomogeneous terms in the equations for the higher modes. 

Derivation of Admittance Relations 

It has been shown (see Eefs. (l2) and (13)) that the solution of Eq. (6) 

can be expressed as the sum of a homogeneous solution A^^^ and a particular 

solution of the inhomogeneous equation A^ ^ as follows: 

P 

Aj,(cp) = KjA^’'^ 9) + A^^\tp) . • ■ (8) 

Using this result a nonlinear admittance relation to be used as a boundary 
condition in nonlinear combustion instability analyses was derived. Noting 
that the velocity potential $ " given by Eq, (5) is a summation of partial 
potentials i ^ where 
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( 9 ) 


■ r =A (<p) 

’ a nozzle admittance relation can te 'written for each of the partial poten- 
tials- This is done by introducing Eq. (8) into Eq. ( 9 ), taking partial 
derivatives with respect to z and t and eliminating K^’ between the 
resulting equations. The resulting admittance relations are given by: 




ih (vtl 
e P J 


r 

P 


= 0 , 


where 



( 11 ) 


( 12 ) 


Equation (lO) is the nonlinear nozzle admittance relation to be used 
as the boundaiy -condition at the nozzle entrance plane in nonlinear stability 
analyses of rocket combustors. The quantities Yp and Fp are, respectively, 
the linear and nonlinear admittance coefficients for the pth mode. The - 
nonlinear admittance, F^, represents the effect of nozzle nonlinearities 
upon the nozzle response, and it is zero when nonlinearities are absent 
(i.e., for the IT mode). 
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It can easily he shown that -when the Mach number at the nozzle entrance 
is small, Eq., (lo) can he e^ipressedj correct to second order, as; 

f 


U 

P 


- Y P 
P P 


- -2 

= - u c r 

e e p 


(13) 


•where and are the cp-dependent amplitudes of 'the ajdLal velocity and 

pressure perturbations respectively. 

in order to use the admittance relation (Eq. (lO) or Eq. (13)) in 

combustion instability analysis, the admittance coefficients Y and T 

P P 

must be determined for the nozzle -under consideration. Phe equations 
governing these quantities are readily derived from Eqs. (6) using the 
definition of F (i.e., Eq. (12) to obtain; 

ir 


H — ^ = 
P dcp 


- M C - M - H C 
P P P P P 


(li^) 


•where 


dF 

H = 

P 








d<? 


(15) 


(l6) 


Calculation of the Mozzle Response 

To obtain the nozzle response for any specific nozzle, Eqs, (l4) and 
(15) are solved in the following manner. As pointed out earlier, the non- 
linear terms vanish for the IT mode (i.e., =0, = 0) and it is only 

necess^ to solve Eq. (l4) to obtain (and hence Y^) at the nozzle 
entrance. Since Eq. (l4) does -not depend on the higher modes, it can be 
solved independently for Xj* Once has been determined both Eqs. (l4) 



and ( 15 ) most iDe solved for. the other inodes. In order to do this, the 
amplitude A^(<p) nmst he determined since Eq. (I 5 ) depends on 
derivatives through Ip(9)‘ Once is hnorvra, A^(t|>) is determined hy 

numerically integrating Eq. (l6) -where the constant of integration is 
determined hy the specified value of the pressure amplitude |p^| (of the' IT 
mode) at the nozzle entrance. The -value of thus found is introduced into 
Eq. ( 15 ) -which is then solved for 

- Since Eqs.(l4) and { 15 ) are first order ordinary differential equations, 

the numerical integration of these equations must start at some initial point 

where the initial conditions are knovm, and terminate at the nozzle entrance 

where the admittance coefficients Y and F are needed. Since the 

P P 

equations are singular at the throat, the integration is initiated at a point 
that is located a short distance upstream of the throat. The needed initial 
conditions are obtained hy expanding the dependent -variables in a Taylor 
series about the throat (tp = O) , 

In Eqs. (l4) and (15) j the quantities H , M , U and I are functions 

p ]p ^ !P 

of the steady-state flcrw variables in the nozzle and these must be coii^iuted 
before perfoming the numerical integration to obtain Cp 3,nd F^. For a 
specified nozzle profile, the steady-state quantities are computed by solving 
the quasi-one-dimens ional isentropic steady-state equations for the nozzle 
flow. Figure 2 shcfws the nozzle profile used in these conputations . All 
of the length variables have been non-domensionalized with respect to the 
radius of the combustion chamber to which the nozzle is attached, and hence 
r =1. At the throat r, , Is fixed by the Mach number at the nozzle 
entrance plane. The nozzle profile is smooth and is completely specified by 
r j +, and 0.., -which are respectively the radius of curva-ture at the 
chamber, radius of cjprvat-ure at the throat and slope of the central conical 
section. The steady-state equations are integrated using equal steps in 
steady-state potential 9 by beginning at the throat and continuing to the 
nozzle entrance where the radius of the wall equals 1. 

A computer program, JTOZAbM, has been developed to numerically solve 
Eqs, (l4) - ( 16 ) and calculate the linear and nonlinear nozzle admittances. . 

A computer code and description of this program is given in Appendix: A. 





COMBUSTION INSTABILITr ANALYSIS 
Combiostioa. Chamber Model 

This section describes the application of the nonlinear nozzle admit- 
tance theory developed in the previous section to the analysis of combustion 
instability in a liquid-propellant rocket' combustor. A cylindrical combustor 
with uniform injection of propellants at one end and a slowly- convergent 
nozzle at the other end was' considered. The liquid propellant rocket motor 
that was analyzed is shown in Figure 3» ihs analysis -of such a motor for a 
linear nozzle response is given in Ref. (ll) . 

The oscillatory flow in the combustion chamber is described by the 
three-dimensional, second-order, potential theory developed in Ref. (H). 

Tn this theory the velocity potential i must satisfy the following nonlinear 
partial differential equation; 


E (§') = $/ + i + i' - ® ' 

" rr r r 2 00 zz tt • 


(17) 


- 2% - V Vt - 

r 


(V-I) ^ H *ee - 


- (Y.1)*; f 


+ vn ^ - t)J - 0 


where Crocco’s time-lag (n — t) model is used to describe the distributed 
unsteady combustion process. In the present analysis the linear nozzle 
boundary condition used in the previous analysis (see Eq. (2) of Ref.‘ ll) was 
replaced by the nonlinear admittance condition given by Eq. (if)). 
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Figure 3* Typical Mathematical Model of a Liguid Rocket Motor. 



Application of Galerkiii Method, 

As Sliming a series espansion of the form, (see Eef . 11 ) : 


= y i = y B (t) cos(m 9 ) J (S r) cosh(ib z) (l8) 

/l_i p 1-1 p m mn p 

p=l p=l 


the Galerkin method "was used tc» obtain approximate solutions to Eq. (l?). 

In Eq. (l 8 ) the radial and azimuthal eigenfunctions are the same as those 

used in the nozzle analysis (see Eq^ 4). Unlike the nozzle analysis where the 

unknown coefficients A (cp) were functions of axial location in the nozzle, 

p 

the unknown coefficients Ep(t) in Eq* (18) are functions of time. The b^ 

appear ing in the axial dependence are the axial acoustic eigenvalues for a 

chamber with a solid wall boundary condition at the injector end and a 

Linear nozzle' admittance condition at the other end. 

The unknown- aapiitudes- • were determined by substituting the 

assumed series expansion (i.e., Eq. (I8)) into the wawe equation (i.e., 

. Eq. (17)) to form the residual E (§ ). Similarly, the series expansion 
» ' c 

was substituted into the nozzle boundary condition (i.e,, Eq. (ID)) to obtain 
the boundaiy residual The .residuals E^(§^) and were 

required to satisfy the following orthogonality condition (see Eef. U) : 




1 

[ E' ($^) Z.(z) @.( 0 )R.(r) rdrdQdz 
c" 3 ■ d d 


(19) 



O ' 


®.( 9 )R.(r) rdrde = 0 
d 0 ' 


j = 1,2, ... IT 



“where the Z. are the • complex conjtigates of the axial acoustic eigenfunctions 
J 

appearing in Eq., (18), and ®. and E. are the azimuthal and radial eigen- 

3 3 

ftmctions respectively. 

Evaluating the spatial integrals inEqs. (19) gave the following system 
of m complex nonlinear equations to he solved for the aii®litude functions, B (t) 


r' 

I 

P=1 


{ 




d^ 


dt 


2 


[ 


CgCdiP) 


- nC. 


,(jjp)] 


dB 
— £ 
dt 


(20) 


+ nC^(d,p) 


d[B^(t-T)] 


dt 




'T* v-> r dB dB 

Z, i iVj.P.liEp - 3 ^ + l>2(3,P,q)B^ 

p=l q=l , ■ 

* * 

* dB ^ dB ^ 

+ l>3(3,P,g)Bj, + D^(j,p,50Bp } = 0 

j = 1,2, N 

In the above equation, the term j ,p)e^p'*^ results from the presence of 
nozzle nonlinearities (i.e» the term involving F in Eq. (lO)). 

The coefficients appearing in Eq. (20) were determined by evaluating 
the various integrals of hyperbolic, trigonometric, and Bessel functions 
that arise from the spatial integrations indicated in the Galerkin oirfcho- 
gonality conditions . These were calculated by the coaiputer. program 
COE5TS3D (Appendix B) . 

The time— dependent behavior of an engine following the introduction of 
a disturbance is determined by speciiying the form of the. initial disturbance 
and then following the subsequent behavior of the individual modes by 
numerically integrating Eqs. (20). Once the time -dependence -of the Individual 
modes is .known, the velocity potential, is calculated from Eq. (l8) . 

The pressure perturbation at any location within the chamber is related to 



"by the following second-order momenttim equa-tion .(see Ref. llj : 


/ 

P 



"z 



( 21 ) 


CTumerical Solution Rrocediire 

Equation (20) is a system of U ordinary differential equations which 
describes the "behavior of the R complex time -dependent functions, 

-B^(t). Beginning with a sinusoidal initial disturhanee, a fourth order 
Runge-Kutta scheme was employed for the numerical integration of this system 
of equations. In the present calculations, a three-mode series expansion 
consisting of the first tangential (IT), second tangential (2T) and first 
radial mode (IR) was used. This is the same series expansion used in the 
stability calculations presented in Refs, (lO) and (H). The numerical 
integration of Eqs. (20) is performed hy the computer program, LCICSD? which is 
described in Appendix C. 

The oscillatory flow in the combustor and nozzle are mutually dependent 
on each other; that is, the combustion chamber analysis requires Knowledge of 
the nozzle admittances, but these nozzle admittances depend on the frequency 
of oscillation and the pressure amplitude, which can only be determined by the 
conibustion chamber analysis. Thus an iterative solution technique is used. 

Ih this procedure, linear nozzle admittances are first calculated for the 
specified nozzle geometry. Rexfc,. the conibustion chamber analysis is 
carried out using these linear nozzle admittances (F = O) , and limit-cycle 
frequency and pressure amplitude of the IT mode at the nozzle entrance 
are determined. This information is then used in the nozzle theory to 
determine the nonlinear nozzle admittances which are xised in the chamber 
analysis to calculate new limit -cycle frequencies and pressure anplitude. 

If the limit-cycle amplitude obtained with the nonlinear nozzle boundary 
condition is significantly different from the limit -cycle amplitude obtained 
with the linear nozzle admittances, new values of the nonlinear admittances 
are calculated and the process is repeated until the change in limit-cycle 
amplitude is sufficiently small. 
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EESULTS AM) Discussion 


AajnxttaKce Coefficients 

Computations of the admittance coefficients hare been performed using 
a three-term series expansion consisting of the first tangential^ second 
tangential and first radial modes. An Adam-Bashforth predictor-corrector 
scheme was used to perform the numerical integration, while the starting values 
needed to apply this method were obtained using a fourth order Eunge-Kutta 
integration scheme. Computations have been performed for several nozzles, at 
different frequencies and pressure amplitudes of the first tangential mode. 

ligure 4 shows the frequency dependence of the linear admittance coeffi- 
cients for the IT, 2T, and IE. modes for a typical nozzle = 20°, r^^ = 1.0, 
r^^ = 0.9234 j M = 0,2). Here, to is the frequency of the IT mode, while the 
frequency of the 2T and JR modes is 2U) due to nonlinear coupling. Hence the 
real parts of the linear admittance coefficients for the 2T and IH modes 
actually attain their peak values at a higher frequency than that for the IT 

mode. The linear admittance coefficients for the IT mode are in complete 

* *1 ^ 

agreement with those calculated previously by Bell and Zinn . 

The frequency dependence of the nonlinear admittance coefficient for 
the 2T mode is shown in Figure 5 with pressure amplitude of the IT mode as a 
parameter. "While the behavior of the linear admittance coefficient depends 
only upon the frequency of oscillations, the behavior of the nonlinear 
admittance coefficient is seen to depend also on the amplitude of the IT mode. 

The absolute values of both T and F. increase with increasing pressure 
anaplitude of the IT mode, which acts as a driving force. It is observed that 
the absolute values of F^ and F^ vary wiuh frequency in a manner similar 
to the absolute values of and Y^. The frequency dependence of the non- 

linear admittance coefficient for the IE mode is shown in Figure 6 with pressure 
aiEplitude of the IT mode as a parameter. 

Figure 7 shows the effect of pressure amplitude upon the magnitude of 
the ratio 'of nonlinear admittance coefficient to the linear admittance coeffi- 
cient for the 2T and IH modes respectively. This ratio, |f/y| , increases 
with increasing pressure amplitude. In the limiting case of |p^| =0, 
the nonlinear admittance coefficient is zero for all frequencies as expected. 
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figure 8 shows the influence of entrance Mach number on the nonlinear 
nozzle admittance coefficients for the 2T and IR modes respectively. Here the 
relative magnitudes of the linear and nonlinear admittances (i.e.j |r/yj) are 
plotted as a function of amplitude of the IT mode. Ei each ease there is a 
significant decrease in jr/Y] with increasing Mach number, thus it appears that 
the in^jortance of nozzle nonlinearities will he smaller at higher Mach nurnbers. 

The effect of nozzle half -angle on jr/Y| for the 2T and IR modes is 
shown in Figure 9* It is readily seen that for 6^ between 15 and 45 degrees 
there is only a slight effect of nozzle half -angle on the relative magnitudes 
of the linear and nonlinear admittances. However, it should, he noted that both 
the linear and nonlinear theories are restricted to slowly convergent nozzles 
(i.e. , small 0^) . 

Figure 10 shows the effect of the nozzle radii of curvature upon the 
quantity |r/Yj for the 2T mode. It is observed that a change in the radius of 
curvature of the nozzle at the throat has an insignificant effect on the 
relative magnitude of the linear and nonlinear admittances. On the other hand, 
a similar change in the radius of curvature of the nozzle at the entrance 
section has considerable effect on the relative magnitude of the linear and 
nonlinear admittances. Similar results were obtained for the IR mode. 

Tn summary, the results obtained in the admittance calculations indicate 
that the magnitude of the nonlinear admittance coefficient is comparable to 
that of the linear admittance coefficient, especially at large pressure an®)!!- 
tudes. To determine if this result has a significant effect upon combustor 
stability, calculations were made for typical liquid rocket coniiustors usxng the 
rinnlnTipa.T admittances. These results were compared with similar calculations 
using linear admittances. The results of this investigation are discussed in the 
remainder of this report. 

Stability Calculations 

Combustion instability calculations have been made using the three mode 
series consisting of the IT, 2T, and IR modes. These calculations have been 
made for different values of the following parameters: (l) time lag t, (2) 
interaction index n, (3) steady state Mach number at the nozzle entrance M^, 
and (4) chamber length-to -diameter ratio L/D. All of the combustors that 
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lia-ve lieen analyzed are attached to nozzles with the foUowxng specifications: 
radins of curvature of nozzle at the combustion chamber, r^^ = 1.0, radius of 
curvature of nozzle at the throat, = 1.0; and nozzle half -angle, 0^-20 . 

Tr> each case, solutions have been obtained with both the linear and nonlinear 

nozzle admittances. 

A typical neutral stability curve is shown in the n-T plane in. figure 11. 
Since it was desired to study the limit-cycle behavior of the motor, the values 
of n a.ni^ T considered were chosen from, the unstable region of this stability 

diagram. 

Limit -cycle amplitudes and wavefoamis were calculated for T = 1.6 
(resonant conditions) for several values of n as shewn in Figure 11. Wall 
pressure waveforms (antinode) are shown for a milddy unstable case (Point A, 
n = 0.52) and a, strongly unstable case (Point B, n = O.70) in Figures 12 and 13* 
Figure l4 shows limit -cycle amplitude as a function of n for t === 1.6. In 
each case both linear and nonlinear nozzle admittances were used in the calcula- 
tions, These results show that the nozzle nonlinearities have only a small 
effect on the limit-cycle amplitude and waveform even for fairly large amplitude 
ins t ab ilit ie s . 

Similar comparisons were made for the off -re sonant values of n and 
T shown in Figure 11 (see points C, D, E, F) . These results alno show very 
little effect of nozzle nonlinearities on the limit-cycle amplitudes for off- 

resonant oscillations as seen in Figure 15 • 

Finally, comparisons of limit-cycle amplitudes are shown for various 
exit Mach numbers in Figure I6 and for various lengfeh-to-diameter ratios in 
Figure IT. Again, limit-cycle amplitudes obtained using the nonlinear nozzle 
boundary condition agree closely with those obtained using the linear nozzle 
boundary condition. 


cowcniD3iir& eemaeos 

A second-order theory and computer program have been developed for cal- 
culating three-dimensional, nonlinear nozzle admittance coefficients to be used 
in the analysis of nonlinear conhustion instability problems. This theory is 
applicable to slowly convergent, supercritical nozzles under isentropic, 
irrotational conditions when the combustion chamber oscillatxons are dominated 
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■by the M mode . Nozzle admittances have "been computed for typical nozzle 
geometries, and results have 'been shc?wn as a function of the frequency and 
amplitude of the IT mode. 

The nonlinear nozzle admittances have "been incorporated into the 
previously developed nonlinear combtistion instability theory, and calculations 
of limit-cycle amplitudes and pressure waveforms have been made to assess the 
importance of the nonlinear contribution to the nozzle admxttanee. These 
results show that nozzle nosalinearities can be safely neglected in nonlinear 
combustion instability calculations if the following conditions are satisfied: 
(l) the amplitude of the oscillations are moderate, (2) the mean flow Mach 
nuniber is small, and (3) the instability is dominated by the first tangential 
mode. Therefore, the linear nozzle boundary condition used in the previous 
nonlinear combustion instability analyses is adequate for most cases involving 
IT mode instability. 
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APPENDIX A 

PEOGEAM NOZADM: A USEE'S MANUAL 

General Description 

Program 3SDZADM calculates Doth the linear and the nonlinear admittance 
coefficients for a specified nozzle. These admittance coefficients are 
required as input for Program C0EFFS3D (see Appendix B) which calculates the 
coefficients of "both the linear and nonlinear terms in the combustor amplitude 
equation (i.e., Eq. (20)). The output of Program NOZADM is either punched . 
onto cards or stored on disk or drum for input to Program C0EPFS3D. 

Program Structure 

A flow chart for Program NOZADM is shown in Fig. (A-1) . The program 
perforins the following operations: (l) reads the input data, (2) calculates the 
steady-state flow quantities in the nozzle, (3) obtains the starting values 
needed to numerically integrate Eqs. (l4) and (l5)j (4) performs the numerical 
integration of Eqs. (l4) and (15) to obtain the desired admittance coefficients, 
and (5) provides the desired output. 

The inputs to the program include parameters describing the nozzle, the 
frequency and pressure amplitude of the fundamental mode, and the various 
control numbers. 

After reading the input, the program obtains the steady-state flow 
quantities at every station in the nozzle by calling the subroutine STEADY. 

This subroutine also calculates the number of station points (NPLAST) in the 
nozzle. 

The evaluation of the admittance coefficients is carried out in stages. 

The work performed in each step depends upon whether or not the nonlinear 
admittances are to be evaluated. If only the linear admittances are required, 
only the equation for C needs to be solved. Thus, the equations govering C 

Jr 

are solved individually for each of the modes in the series expansion. On the 
other hand, if the nonlinear admittances are also required the equations 
governing the linear admittance for the fundamental mode (£ ^^) and the amplitude 
of the fundamental mode (A„) are first solved to obtain these quantities at 
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I'igure A-1. Elow Chart. 












every station in the nozzle. In the siibsequent steps, the equations for C 
and r for each of the remaining modes are solved. 

Input Data 

A precise definition of the input data required to run the computer 
program is given helow. The input is given through three data cards. 3h the 
description of the cards helow, the location number refers to the columns of 
the card. "l” indicates integers and "f" indicates real 3iumbers with a 
decimal point. For the I formats, the values are placed in fields of five 
locations while a field of ten locations is used with the "f" formats. In 
either case, the numbers must he placed in the rightmost locations of the 
allocated field. 

No. of 


Cards 

Location 

Type 

Input Item 

Comments 

1 

1-10 

F 

CM 

Mach number at the nozzle 
entrance 


11-20 

F 

ANGLE 

Nozzle half -angle 


21-30 

F 

ECC 

Radius of curvature of the 
nozzle at the entrance 


31“40 

F 

RCT 

Radius of curvature of the 
nozzle at the throat 


la-50 

F 

GAM 

Ratio of specific heats 

1 

1-5 

I 

NOZNll 

If 0: nonlinear'’ admittances 
are not evaluated 

If 1: nonlinear admittances 
are evaluated 


6-10 

I 

NOUT 

Determines output 


If 0: only printed output 

If 1: printed- and stored on 
disk or drum (output device 
number 7) 

If 2t printed and cards 
punched in a format suitable 
for the program C0EFFS3D 



No of 
Cards 


1 


Location 

11-15 

I 

Input Item 

Comments 

lEXTN 

If 0: no extension section 
If L: an extension section 
is present. 

16-25 

F 

EXTNSN 

length of the extension 
section; canit if lEXTN - 0 

1-10 

F 

WC 

Frequency of oscillation 

11-20 

F 

PIAMPL 

Pressure anplitude of the 
fundamental mode. Omit if 
only linear admittances 
are needed. 


The nozzle parameters ANGLE j HCC and HOT correspond to 6^5 r^^ and r^^ 
in Nig. 2. For lEXTN = 1, the integration of Eqs. (l4) aJ^d (I5) is continued 
heyond the nozzle entrance plane to a length EXTNSN within the combustion chamher 
When NOUT = 1, the -values of the necessary aamittance coefficients are stored 
on disk or drum (device number 7) ia a format suitable for nuiput to program 
C0EFFS3D. If, instead of providing this data to program C0EFFS3D through data 
file 7, it is desirable to provide punched cards only, NOTjT should be 2. 

Again the format is such that these cards can be fed to program C0EFFS3D 

directly. , 


Steady-State Quantities 

The subroutin^^TEADY is called to evaluate the steady-state quantities 
in the nozzle. This subroutine first calculates the radius of the nozzle at 
the throat necessary to obtain the specified Mach nurber at the nozzle entrance. 
The steady-state flow quantities at the throat are determined by the choking 
conditions. Starting with these values, .the steady-state flow quantises at the 
other stations in the nozzle are calculated by humerically integrating the 
steady-state equations starting from the throat. The subroutine EKSTDY deter- 
mines the values of the steady-state velocity near the throat using the 
Runge-Kutta scheme. These values are needed to start the Adam’s predictor- ^ 
corrector scheme for integrating the steady-state flow equation. The numerical 
integration is performed by the subroutine UAIAI©. Starting slightly upstream. 
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of the throat, the numerical integration is continued till the nozzle entrance 
is reached (radius of the nozzle R = l) . The arrays V and C contain the ' 
steady “-state velocity and speed of sound respectively. 


Coefficients 

The complex coefficients that appear in the nozzle admittance equations 
are evaluated in the program hy calling the subroutine COEFFS. These coeffi- 
cients contain certain integrals involving trigonometric and Bessel functions. 
The subroutine HJTGRL sets up arrays for these integrals. 


integrals 

The necessary trigonometric integrals are determined by the subroutine 
IKTGEL itself. Denoting 

® (0) = cos(m 0), 

p P 


the integrals are as follows: 
ALPHA (1, p) 

Alim (2, p) 

ALPHA (3, p) 

ALPHA (4, p) 

ALPHA (5, p) 


2tt j- -.2 

J [®pCe)J ®i(e) 


o 

2rr 


J [®'(9)] ®i(e) aa " 


O 

2tt 


J ®"(0) ® (0) ® (9) 


o . 
2tt 


P' ' p' ' r 

.2 






d0 


o 

2tt 


f ®"(e) ® .(9) de 

« T) T) 
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The integrals involving Bessel functions are as follows: 


Here 


ju 2 

BETA (1, p) = ^ '^r 

o 

1 2 

BETA (2, p) = J f ^ 


1 


BETA 


(3, P) = J ^ ^ 


BETA 


(4, p) = J Rp(r) K^(r) E^(r) r 
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BETA 


(5, p) = J R^(^) K (^) ^ 

^ P 


BETA 


2 

(6, p) = j ^ 


BETA 


JU 

(7, p) = J RpW Rp(^) ^ 


BETA 


(8, p) = J Kp(r) Rp(r) r dr 


1 


BETA 


.2 1 ^ 
— dr 


r r A “ 
(9, p) = J L^p(^)j ^ 


o 


R fr') = J IS rl where m and n are the transverse mode numbers 
p^ -mL mn J 



These integrals of Bessel functions are obtained from the functions 
RADI and EAB2. BM>2 provides the first five integrals while EABl provides 
the last four integrals, Simpson’s integration scheme -is .used in these 
function subprograms to evaluate these integrals. The values of the Bessel 
functions of the first kind are obtained using the subroutine JBES (see Ref. 17 

Integration of the Differential Equations 

For the numerical integration of the differential equations, a fourth- 
order Adam-Bashforth predictor-corrector scheme is employed. The necessary 
initial values are obtained by using a fourth-order Eunge-Kutta scheme near 
the throat. The Runge-Kutta integration is performed by subroutine RECTZ. 

The predictor-corrector integration is performed by subroutines TABAmS and 
ZADANB. The values of the dependent variables are stored in the array Y and 
their derivatives are stored in the array DY. The integration is continued in 
steps of BP in the axial variable (steady-state velocity potential) till the 
combiAstion chamber is reached. 

After the numerical integration of all the differential equations is 
completed, the admittance coefficients are evaluated. .AMPL (j) and 
PHASE(j) are the amplitude and phase of the linear admittance coefficient for 
mode J. GH'OZ(J') is the complex, nonlinear admittance coefficient for mode J. 

Output 

The output of the program NOZADM contains ttTO sections. 

In Section 1, the parameters of the nozzle being analyzed are printed 
out. The output of this section occupies only one page and is essentially a 
print out of the input data. The parameters, which are printed are: 
the Mach number at the nozzle entrance (CM), the specific heat ratio (gAM), 
the nozzle half -angle (ARGIE), the length of the extension section, if any 
(EXTMSN) , the radius of curvature of the nozzle at the throat (ECT) , the 
radius of curvature of the nozzle at. the entrance (rcc) , and the number of 
stations in the nozzle (hPIAST) . Section 1 is printed for any value of the 
control number IfOBT. 
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Section 2 contains the nozzle admittance coefficients. Depending on 
the value of the control number STOUT, Section 2 is printed, stored on disk or 
drum or punched onto cards . These three modes of output will now he discussed 
individually . 

Printed output ; The control number NOUT for this mode is 0. The 
printed output appears on one page and contains hoth the linear and nonlinear 
admittance coefficients. For each coefficient, the real and imaginary parts as 
well as the magnitude and phase are printed out. If nonlinear admittance 
coefficients are not calculated hy the program (iTOZHLl = O) , zeros are entered 
in the spaces for the nonlinear coefficients. 

This mode of output is inconvenient to use for instability analysis 
since it would then he necessary to manually punch all the input cards for the 
program C0EFFS3U . 

Disk or Dmrn Storage : The control number NOUT for this mode is 1. 

When disk or drum storage (like the FASTKADD System on the UNIVAG 1108) is 
available, this is the most convenient means of storing the output of 
Section 2. The necessary admittance coefficients are stored in a format suit- 
able for input to the program C0EFFS3D. The device number for this output 
is 7. The control statement needed to request the disk or drum storage on the 
computer depends on the computer facilities being used. 

Punched Cards ; WOUT for this mode is 2. This mode of output is the 
simplest way to run the instability program. The cards containing the 
necessary admittance coefficients are punched by the computer in a format 
suitable for use with program C0EFFS3D, which is the next program to be executed. 
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FORTRAN Listing 
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PROGRAM NOZAEM ♦♦♦sSi*************************** 

THIS PROGRAM EUAL OATES THE LINEAR AND NONLINEAR ADMITTANCES 
OF A SPECIFIED NOZZLE. 

THE FOLLOWING INPUTS ARE 'RESUIRED : 

CM IS THE MACH NUMBER AT THE NOZZLE ENTRANCE. 

ANGLE IS THE SLOPE OF THE MIDDLE SECTION GF THE NOZZLE. 

RCC IS THE RADIUS OF CURVATURE OF THE NOZZLE AT THE ENTRANCE. 
HCT IS THE RADIUS OF CURVATURE AT THE THROAT- 
GAM IS THE SPECIFIC HEATS RATIO. 

NOZNLl DETERMINES WHETHER THE NONLINEAR ADMITTANCES ARE TO 

BE evaluated: 

NOZNH » 0 NOT EVALUATED. 

NOZNLl = 1 EVALUATED. 

ROUT DETERMINES THE OUTPUT: , 

NOUT = 0 PRINTED OUTPUT ONLY- 

NOUT * 1 PRINTED AND WHITTEN INTO A FASTRAND FILE. 

NO.UT * e PRINTED. AND ADMITTANCES PUNCHED INTO CA.RDS. 

I EX IN' DETERMINES IF THERE IS AN EXTENSION SECTION 
lEXTN *=0 NO EXIENSIOK SECTION. 

lEXTN = 1 THERE IS AN EXTENSION SECTION. 

EXTNSN IS THE LENGTH OF THE EXTENSION SECTION. 

WC IS THE FRE&UENCY OF THE FUNDAMENTAL MODE. 

PIAMPL IS THE PRESSURE AMPLITUDE OF THE FUNDAMENTAL MODE. 


COMMON /X 1 / CM* ANGLE* RCC* RCT* G AM* G* RT* DP 

1 /X2/T*H1*RS*NPLAST*NEND*IEX‘1N - 

/Xa/WC* SVN*IP*MODE*NU*KPC35 
/XA/RUC7>*RDUt7)*ZIHRl*GTHRl 
/X5/U< 1000>* DU< 1000)* C( 1000)*RVC 1000) 

/XG/AFK* AFN1*AF-N2 
/X7/ALPHAC5*3)* BETA(9*3) 

✓X8/ZRKC 1000) 

AFN( 1000)* AFNK 1000)* AFN2C 1000)* ACHKBH* CONST* 
CC< £5) * CCl C25)* CFh* CFM* CFN* CGRP J * CGRP2* 

INHKG* 1NHMG1*ZTHR*ZTHRI*AH* AHl*OTHR* 61HB1* 
ZETA* TAU*LINAEM*ZRK*GN0Z<3) 

DIMENSION 6CA)*CP(4)*YCA)* DY( A* A)* SMNC 3)*1STEPC3)* 

I NAME<3>* PHASEC3)* AMFLC3) 

DATA <NAME<MOr>E)* MODE - 1*3) /8H1 T* 2H2T* 2H IR/ 

1 CSMNCMODE)* MODE « i.3> / 1 .8A1 18* 3.05A2A* 3.83 J7 IZ 


READ <5*5005) CM* ANGLE* BCC*RCT*GAM 
READ <5*5010) NOZNLl* NOUT* lEXIN* EXTNSN 
READ <5*5015) WC* PJA!«1PL 
'GKINl a GAM - 1. 

GFLl a gam +1. 

DP »= -0.002 

I STEP « 1 : INTEGRATE FOE ZETA ONLY* 



48 



c 

c 


10 


15 


c 

c 

c 

c 


c 


c 

c 

c 

c 


20 

S5 

C 


C 


C 


ISTEP « £ t IRTEOPATE. POP ZETA i AH. 

ISTEP t: 3 9 INTEGHATE FOK ZETA * GAKWA. 

IF <R0ZML1 .E0« 1) GO TO 10 

1STEF<I) = 1 

ISTEFCS) r= 1 

ISTEP<3> =1 

GO TO 15 

ISTEFCl) B 2 

lSTEF<e> * 3 

ISTEP<3> B 3 

CONTINUE 

KP< 1) «= 1 

KF<2) t= 2 

KP<3) e 2 ■ 

OBTAIN STEACY” STATE OUANTITIES IN THE NOZZLE . 
CALL STEADY 


PRINT OUT THE -NOZZLE PARAKETEPS. 
UBl TE <6> 1005) 

VRITE <6* 1010) CM 
VRITE <t* 1015) gam 
VEITE C6*1020) ANGLE 
VRITE <6^ 1025) EXTNSN 
VRITE <6# 1030) RCT 
VRITE (6.. 1035) RCC 
VRITE (OjIOAO) WFLAST 


NEND * NPLAST 

IF CIEXTN .NE* 1) 60 TO 25 


DETERMINE NUMBER OF STATIONS IN THE EXTENSION REGION^ AND 
DEFINE STEADY- STATE QUANTITIES IN THAT REGION. 


UEXT = U(NPLAST) 

NEND = NPLAST - « EXTNSN * UEXT ♦♦ *5) / DP 
DO 20 NP = NPLAST/ NEND 
UCNP) B U< NPLAST) 

C(NP) * CCNPLAST) 

DU<NP) B DU< NPLAST) 

RW(NP) « RV< NPLAST) 

CONTINUE 

CONTINUE ORIGINAL PAGE IS 

IF (NEND .GT. 1000) GO TO 550 OP POOR QUALET^ 

CALL INTGRL 
£RTR«(RT*RCT>**. 5 


ACHMRR «= CMPLX (PIAMPL 
IF (NGUT .EQ. 0) WRITE 
IF (NOUT .EO. 0) WRITE 


/ <WC*6AM> /O.) 
(6*1050) VC*P1AMFL 
(6* 1055) 


DO 500 M0DE=1*3 

IPbISTEP(MODE) 

SVN«SKN(«ODE) 
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SVKH=SVN/BT 


C 

VALUES SECTI0N+*++>t'*+=i'*******+1=+***+++it********>»* 
C 

p=0. 

ABB •: 1* 

AHI c 0* 

AH = CMFLX CAHR^AH1> 

UP 5= Uf I ) 

CP >= C( l> 

PUP = DUU> 

RVF •= RWCi> 

CALL COEFFS <UPi DUP^ CF*BVP^ CC) 

CFH *= CCCi) 

CFM = CCC85 + CC<6) 

CFN = CCC3) + CC(4) + CC<53 + CC<7) + CC<8> 

peri VATIVES OF THE COEFFICIENTS AT THE THROAT^*#’*'**** 

EVALUATE DERIVATIVES OF LINEAR COEFFICIENTS* 

XR « - A*/'C6PLI ♦ SRTR> 

CFHl = CKPLX CXR>0*> 

XR i: - (fiA. ♦ A* * GAK> ✓ <CFL! ♦ 3* ♦ RT * RCT> 

XI c - g, + VjC * KPCM0DE5 / CGFLl * SRTR) 

CFMI = CMFLX <XR^XI> 

XR = ♦ S.*6N1N1 + CBETA <8*K0DE) + BETA' CT^NODE) + BETA C9»N0DE) 

I * ALPHA <5* NODES / AJ-FHA MOPES S / <GFL1 * HT ♦ RT 

8 + SRTE BETA C6>K0PF)S ' 

XI «= 12 + 8*GAM) ♦ WC * KP(MODE) ♦ GMINl / <3.*GPL1 * RT*BCT) 

CFNl *= CMPLX CXR^XIS 

SET UP VALUES AT THE THROAT ■ BT, TAYLORS EXPANSION 

STARTING VALUES' FOR ZETA 
ZTHR *: - CFN / CFM 

ZTHRl e - (CFMl ¥ ZTHR + CFHI * ZTHR * ZTHR + CFNIS / <CFHl + CFMS 
ZRKCl) s ZTHR 
C 

IF (M0P£*NE«1> GO TC 110 
AFKC 1 ) c AH 

AFNKl) « AFNC13 =9= ZTHR 

AFN8(1S = AFNKIS ♦ ZTHR + AFN<U ♦ ZIHRi 
110 CONTINUE 

GClS = HEAL CZTHRS 
GC8S fc AIMAG CZTHRS 
DY Cl»l) e real CZTHRIS 
DY C£j 1) *= AIMAC CZTHRIS 
CO TO < 120* 130# 1 AOS ^ IP 
130 GC35 = AHR 
GCAS *: AHI 
AHI c AH * ZTHR 
DY C3*1S = REAL CAHIS 
PYCA#1S = AIMAG CAHl) 

GO TO 180 
I AO CONTINUE 
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CGRPI *= 
CGHP8 » 
INHMG ■ 

S 


CCC13) + CCCiiil + CCC195 + CCC233 + CC<8^) ♦ CC<25) 
CCC103 + CCCin + CC<173 + CCC203 + CCCS13 ♦ CCt82J 
-CC<18> * * AFN2C1) - CC«18> ♦ AFNU13 * AFNSCU 

-CCCC93 + CC<15)3 * AFN1C15 ♦ AFNICI) - CCRPl AF«< P * 
AFNJCP “ C6BP2 * AFNCIJ * AFW< P 
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EVALUATE DEPIUATIVES OF NON-LINEAH COEFFICIENTS* 

AlBl e ALPHAt I>MOOE) * BETA.< PWODE) 

A£BS c ALPHA<£>l»)ODE> # BETA<S*MOrE> 

AXB3 « ALFHA< UMODE) ♦ BETA<3>MDCE) 

AAB6 e ALPHAS AjNOBE) * BETA<6>«0DE> 
to £6 0 » 1/25 

CCKO) e CKPLX <0./0o) 

XR ® - <2**A1BI JOt WC> / <AAB6 ♦ GFLl * SBTR> 

XI * XB 

CCK9 3 *= CKPLX (XR,XI) ^ 

XB * * CA* * AlBl) / C 3. 1415927 ♦ GPLl * SBTR * A4B6) 

XI * - 3CB 

CCl <12> » CMPLX <XE/XI) 

XR e - A1B3 / <GFLl * RT * RT + SRTH ♦ A4B6> 

XI « - XR 

CCl <133 = CMPLX CXR/XI3 

XR » - A2B2 / <GPL1 * BT ♦ RT * A4B6 * SRTR> 


XI = - XR 

CCl <U> “ CMPLX CXR/XI3 
XR * “ AlBl * C3e*GFLl ♦ SRTR + GMINl 
I <£- * RT « RCT * GFLl GFLl # 

XI « - XR 

CCl <153 « CMPLX CXB/XI3 
XB “ A1B3 » <9* “ 2*»GAM ” GAM*GAK3 
I ♦ A4B63 


♦ <18 
A4B63 


.+GAM>> / 


/ <!£• * RT**3 * RCT >1' GFLl 


XI e - XR 

CCl <163 = CMFLX <XR/XI3 
XR = A£B£ ♦ <9* - £**GAK 
1 A4B63 


GAM*GAM) / <!£• ♦ RT*«3 + RCT * GFLl 


XI » - XR 

CCl <173 •= CMPLX (XR/X13 
XR *= - <GM1NI * VC # AlBl) 
XI « XR 

CCl <183 >= CMFLX <XR/XI> 

XB • <GMIKl * <6.+6AM) # 

I 3Sc AAB63 


/ (GFLl »)' SRTR * A4B63 


VC V AlBl) / <3* + GFLl ♦ BT ♦ RCT 


XI = XR 

CCl <193 = CMPLX CXB»XI) ^ 

XR B - <GMINl * ALPHA < 1/MODE) * (BETA <4/MODE> • 
1 / <GFL1 * BT >»= RT * SRTR * A4B63 

XI *= •> XR 


CCl <233 = CMFLX CXRrXI) 

XB » - <GMINl + ALPHA < l/MODE) * BETA <5*M0DE> * 
/ (GPLl * RT ♦ BT * SRTR ♦ A4B6> 

XI *= - XR 

CCl <243 = 03PLX CXB/XI3 wr.r,trx-. 

XR «s - <6MINl ALPHA < 3/MODE) * BETA <2/M0DE>) 

[ / CGPLl # RT RT * SRTR * A4B63 

XI « - XR 


EETA<5/MODE3)3 


£•> 
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CCl <85) e CMPLX CXR»XIJ 
C 

INHMGl c - AFNEt n * AFN£< U * CCC123 - AFNIC15 * AFN2C 1> * 

1 CCC<18> ♦ CCUlE) + a-*f'CC<9) + E.*CC<15)) - AFN2(i) 

8 ♦ AFK(U * CCC1U8) + CGFiFn - AFNKJ> * AFNICU + 

3 CCCK9) + CCKIS) + CGBPl) - AFKICI) ♦ AFM( t> ♦ 

A CCCK13) + CCKIA) + CC1C19> + CC1CS3) ♦ CCJC2AJ 

5 + CCIC25> + e.*CGRPE> * AFN( 1) * AFN< U * <CCICJOJ 

6 + CClUi> + CCK17> + CCHEO) ♦ CCHED ♦ CCICSEW 

C ■ ■ 

C STARTING VALUES FOR GAMMA 

GTHR e - INHMG / < CF * CFM) 

GTHEl e <-CF * GTHR * CCFHl ♦ ZTHR + CFKl) ♦ CGMINl ♦ *5 » A* / 

1 < GFLl + SBTR>i * GTHR * CCFHl + CF«> “ IKH«Gl> / 

S < CP * CFHl + CF ♦ CFM) 

C 

GC33 = REAL CGTHR) 

6<A> = AIMA6 CGTHR) 

EY C3#l) = REAL CGTHRl) 

DY <A>1) « AIKAG CGTHRl) 

120 CONTINUE 
C 

c ♦♦♦♦♦♦♦♦♦♦♦♦♦NUKERI cal COMFUTATI OUS******^*********’t***’^*********^ 

c 

c hunge-kotta integration to provide initial values 
C FOE predictor- CO ERECTOR IKTEGFJiTION 

C 

DO 30 IRK * 2» A 

CALL RKT2<DP^P#6#GP»1RK) 

p»p+'dp 

ZR«6( 1) 

ZIt=G(8) 

ZBKCIBK) = CMFLX CZR«ZI) 

DY< 1/IRK>bGP< U 
DY<e^IRK>*=GP<8) 

GO TO C 150> 160> 170)> IP 
160 AHR « 6C3) 

AHI » GCA) 

DYCS^ IRK>«=GPC3) 

DYCA>IRK)s*GP€A) 

IF CMODE*NE.l) GO TO 162 
AFK URK) « CKPLX CGC3)*G<A)) 

AFNl CIRK) = CMPLX C GPC 3>* GPC A) > 

AR2 =» 6(t>*6P<3) - GtS)*GP<A) + GPC1) + G<3) - GPC8)*GCA) 

AI2 = GC2>*GPC3) + GC1)*GPCA) + GFCE)*6(3) + GFCD + GCA) 

AF NSC IRK) e CMPLX (ARE^AIE) 

168 60 TO 150 

170 CONTINUE 

GAME = 6C3) 

GAMJ c GCA) 

DYC3>IRK) « GFC3) 

DYCAjIFK) «= GP<A> 

150 CONTINUE 
30 CONTINUE 
YC1)«ZR 
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YC2>e21I 

GO TO < 180# 190*e00>* IP 
190 y<3) « AHH 

Y(/j) « AHI 
CO TO 180 
200 CONTINUE 

YO> Gi^R 
Y<4) • Gmi 
180 CONTINUE 
C 

C PREDICTOR- CORRECTOR INTEGKPiTtON 

CALL E ADAMS ( DP# P# Y# DY# I TOBZ > 

C 

C 

C CALCULATE LINEAR A£«ITTANCE COEFFICIENTS. 

UE B U<NEND) 

CE e C<NEND> 

RHOE «= CE ♦* Cl./GK1N1> 

FR ■ VfC * KFCMODE) 
r ® UE ♦* .5 /■ <FB'^GAM) 

IF (I TORE .EQ. U GO TO 3S 
ER«YCl> 

EIcY<8> 

ZETA « CKPLX CZR>ZI) 

LINAD«1 = F + CMFLXCO.#!*) * ZETA 
GO TO AO 
35 TRc YC1> 

TI YC2> 

TAU * CMPLX (TR#TI> 

LlNAtM B F * CKPLXCO.#l.> / TAU 
AO CONTINUE 

YR «= REAL (LINAD«> 

YI = AIMAG CLINADM) 

YMA6 » CABS CLINAEM> 

YPHASE = ATAN8 CY1#YR> * 180. / 3.1A1S927 
AMFLCMODE) = YMAG 
PHASECMODE) « YPHASE 
C 

GO TO Cai0#220#230)# IP 
220 AHR = YC3> 

AHI e Y(A> 

IF CMODE .NE« 1) GO TO 210 
CONST *= ACHMBR / AFNCNENDl 
DO 50 NP *= 1#NEND 

AFNCNP) = CONST * AFNCNP3 
AFNICNP? == CONST * AFNKNP) 

AFN2CNP> « CONST ♦ AFN2CNFJ 
50 CONTINUE 

C 

C NONLINEAR AEMITTANCE COEFFICIENT IS ZERO FOR IT MODE* 

GAMR B 0. 

GAMl «= 0. 

GMAG * 0. 

GFHASE B 0* 

GBYY e 0.0 



GNOZCU s (0«0>0>0) 

C 

GO TO SIO 
230 CONTINUE 
C 

C CALCULATE NONLINEAR AW1ITTANCE COEFFICIENTS. 

GAMS = V<3> 

GAMI e Y<4) 

GMAG e <GAMK ♦ GAMR ♦ GAMI ♦ GAMI) ** . 5 
GPHASE * ATAN2 CGAMI^GAMR) * 180. / 3.1415927 
GBYY S' CAES (CMFLX <CAKR#6AMn / LINADM) 

GNOZ<MOt£> = CMPLX< GAME* GAMI > 

C 

210 CONTINUE 

IF CNOUT .EQ. 0> WRITE C6> 10605 NAME(MODE>^ YE# YI# 

I YMAG# Y PHASE# GAMR# GAMI# GMAG# GPHASE# GEYY 
SOO CONTINUE ' 

510 CONTINUE 
520 CONTINUE 
550 CONTINUE 

IF CNOUT .EG. 0> 60 TO 560 
to 570 U »= I# 3 

IF CNOUT .EG. 15 WRITE C7# 70055 J# AKFLCJ5# FHASECJ) 

IF CNOUT .EG. 25 PUNCH 7005 U# AMFLCU)# PHASECJ) 

570 CONTINUE 

IF CNOZNLl *EQ. 05 GO TO 560 
DO 580 J »= 1# 3 

IF CNOUT .EG. 15 WRITE (7#70055 J# 6N0ZC0) 

IF CNOUT »EQ. 2) PUNCH 7005 d# GN0ZtJ5 
580 CONTINUE 
560 WRITE <6# 10655 
C 

c ***#♦**♦♦* HMfe****** SEAD format SPECIFICATIONS *♦*♦*♦#»*♦*♦♦♦««***♦ 
C 

5005 FORMAT C6F10.0) 

5010 FOPi^AT <31 5# FI 0-05 
5015 FORMAT C2F10.05 
C 

c 

c WRITE FORMAT SPECIFICATIONS *♦**♦♦*♦♦*<■♦*****+♦* 

C 

1005 FORMAT < IH 1 #////✓✓///# ASX# 1 7H****************** /# 45X> 

I 17KN0ZZLE PARAMETERS# /# ASX# I THtX'************!*:**# //////✓) 

1010 FORMAT C 1H0#25X#”KACH NUMBER "/FA. 25 
1015 FORMAT C IHO# 25X# “GAMMA = •*#FA.2> 

1020 FORMAT C IHO# 25X# “NOZZLE ANGLE *= “#F5«25 

'1025 FORMAT C IHO# 25X# “LENGTH OF EXTENSION SECTION *= “#FA.2> 

1030 FORMAT < IHO# 25X# “RADI US OF CURVATURE AT THE THROAT *= “#F7.55 

1035 FORMAT C IHO# 25X# “RADI US OF CURVATURE AT THE NOZZLE ENTRANCE .= *'# 

I F7.55 

lOAO FOF^AT C 1H0#25X#“NI»1BER OF STATIONS IN WE NOZZLE b “>I45 
1050 FORMAT C IH 1* ////# A6X# 18H******#***»******** /> A6X# 

1 16KN0ZZLE Ami TTANCES#/# A6X# l5H»****’i'+>i'*«**+****** 

2 20X#"FREGUENCY = “# F8. 6# AOX# "PRESSURE AMPLITUDE = “#F6«A5 
1055 FORMAT <///////# 5X# “MODE"# | OX# 2RYR#9X# 2HYI #9X# “YMAG">9X# “YFHASE"# 
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OEIGINAL PAGE IB 
OP POOR QUALITY 



o o o 


1060 

1065 

7005 


1 U34>2HCfU9X#£HGI*9X#4H6«Ae# 10X#6H0PHASE» 1 3X/ 3H6/Y# //) 

FOEMAT < lH0#5X/A8»eX>3FlS*4*F16»A#3F12.4»eF16.4> 

FOEWAT ClHl) 

FO»0AT CIS>aF10.5> 


♦♦♦♦**♦**♦♦♦**•*»*♦*♦**♦♦♦♦♦**♦♦**♦♦♦***♦♦♦♦♦*♦♦♦♦ 


STOP 

END 
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SUBROUTINE STEADY 


C 

C THIS SUBROUTINE EVALUATES STEADY-STATE QUANTITIES IN THE NOZZLE# 

C 

C NOZZLE PROFILE AND FLOW PAR/«<ETERS ARE PASSED TO ‘THE SUBROUTINE 
C THROUGH THE COMMON BLOCKS XI AND X2» 

C THE SUBPROGRAM PROVIDES THE OUTPUT THROUGH COMMON BLOCK X5* 

C U IS THE SQUARE OF THE STEADY-STATE VELOCITY! 

C DU IS THE DERIVATIVE OF U WITH RESPECT TO STEADY-STATE POTENTIAL i 

C CIS THE SQUARE OF THE SPEED OF SOUND! 

C BW IS THE RADIUS OF THE NOZZLE. 

C THESE OUTPUT QUANTITIES ARE STORED IN THE RESPECTIVE ARRAYS AT 

C INTERVALS OF DP IN P « STEADY-STATE POTENTIAL). 

C 

C 

COMMON /XI/ CM*ANGLE.RCC> RCT.GAM. Q.RT.DP 
COMMON /XS/ T#R1^R2.NFLAST*NEND> I EXTN 
COMMON /XA/ RU€7)>BDUC7).ZTHR1^ GTHBl 
COMMON /X5/ UC 1000)*DUC1000)j C( IOOO)» RW< 1000) 

C 

T*= 3* 141 sgR?* ANGLE/ 180* 

RT « <CM**.5) ♦ C<l.+CGAM-l.)*CM+*2/2.) ** (C-C^W-1«)/ 

1 (4»<GAM“l))))*(t2/<GAM+l)> C < -G/M- I ) /C /!•* ( GAM- 1 ) ) ) ) 

SRTR = CRT*RCT) *«*5 

Q *s C.25+RT) * < C2./<6AM+ !•) ) CtCAM+l*) / < GAM- t ) ) ) ) 

B1 «= RT+RCT*< 1*-C0S(T>> 

R2 » l.-RCC ♦ <I«-COSCT)) 

R=RT 
P= 0. 

BW( 1) = BT 

UO) *= 2./CGAM-*-!.) 

Ru< n uc 1 ) 

C< 1) = UCl) 

DU<1) •= 4./<(6AM+l.)*SRTR) 

RDUU) « DUtl) 

G s UU> 

DO 30 1<»2.7 

CALL RKSTDY <P,G.GF) 

P *= P + DP/2* 

RU(I) “ G 
RDUa )*=6P 

IF <I .EG. 2*<I/2>> GO TO 30 
NF = (l+D/2 
U(NP) a RUCl) 

DU(NP) = RDU<I> 

C(NP) ** l.-<GAM-l)-s=U<NP)*»S 

RW<NP> = Q**C<CCNF>) +♦ 1 . /( 2.*< GAM- 1 . ) )) ) 

I *CU<NP)**-.25)*A. 

30 CONTINUE 

CALL UADAMS (P) 

RETURN 

END 
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SUBROUTINE RKSTUYCP>6#DtJMJ 


THIS SUBROUTINE PERFORMS A FOURTH ORDER RUNGE-KUTTA INTEGRATION 
TO OBTAIN STARTING VALUES OF STEABY-STATE VELOCITY FOR THE 
PREDICTOR-CORRECTOR METHOD. 

PIS THE CURRENT VALUE OF THE STEADY-STATE POTENTIALS INPUT* 

G IS THE SQUARE THE STEADY- STATE VELOCITYS INPUT AND OUTPUT. 

AS OUTPUT/ G Is THE VALUE AT THE NEXT STEP* 

DUM IS DERIVATIVE OF THE SQUARE OF STEADY- STATE VELOCITYS OUTPUT- 
BUM IS OBTAINED BY CALLING SUBROUTINE EKUBIF* 


COMMON /XI/ CM/ANGLE/BCC/RCT/GAM.Q/RT/BP 
DIMENSION A<A)/FE<A> 

C 

A(l> e 0* 

A(R> B 0*5 
A(3> B 0*5 
A( A> B 1. 

H B BP/S. 

FR*P 

GR«*G 

CALL RKUDl F< PR/ GE/ DUM > 

FZC.n c D 01 

BO 30 Ib2/A 

PR e P+A<1>*H 

GR c G+'ACl )>»:H»FZCI-1) 

CALL RKUDl F CPR/GB/BUM> 

FZCn = DUM 
30 CONTINUE 

G B 6 + H* CFZCl) + 2+CFZ<2>+FZ<3)> + FZ<A>J/6« 

CALL RKUDl FC PE/ G/ DUM 3 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

10 

15 


22 

25 

30 

35 

40 

45 


50 

20 


SUBROUTINE RKUDIFt GF> 

THIS SUBROUTINE EV/U.UOTES THE DIFFERENTIAL ELEMENT IN THE 
HUNGE-KUTTA INTEGRATION SCHEME FOR SOLVING THE FGUATION FOR SQUARE 
OF STEADY-STATE VELOCITY. 

PIS THE VALUE OF STEADY- STATE POTENTIAL AT THE STATION# 

WHERE DIFFERENTIAL ELEMENT IS SOUGHT! INPUT. 

GIS THE VALUE OF THE FUNCTION AT F5 INPUT. 

GP IS THE REQUIRED DIFFERENTIAL ELEMENT. 

COMMON /XI/ CM#ANGLE#RCC#RCT#GAM#e#RT#DP 
COMMON /X2/- T#Rl>R2#NFLAST#NEND#IEXTN 
COMMON /X3/ VC#SUNj1P#K0DE>NU#KF<3) 

IF CF5 IS# 10# IS 

GF e 4./ CCGAM+1.3 ♦ t(RCT*RT) **.5)) 

GO TO 20 

C *= l-tGAM-l. >+G+*5 

B e e*<<C> ** <-l-/(E.*<eAM-1.5?5 J * <G*>t=-.8£> * 4* 

IF CR-1.) 22#22#50 
IF «R-R1) 25# 30# 30 

DR >: -<<2.*RCT*(R-RTJ - <R-RT) ♦ <R-RT))**.5> / CRT+RCT-FO 
GO TO 45 

IF <E-R2) 35# 40# 40 

DR -TANtT) 

GO TO 45 

DR «= C<2.*BCC*( 1-R5 - <B-1>*CR-15? **.5) / Cl.-R-RCC) 

DU »» -<G**.75)^CC**<<2.*GAM-l) / <2.*< GAM- 1. > >Y > >' 

1 (Q*C I.-CGAK+l. > ♦ G+-53) 

GF e DU*ER 
GO TO 20 
GP *= 0. 

RETURN 

END 
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SUBROUTINE UADAMSCPJ 


C 

C THIS SUBROUTINE CARRIES OUT A MODIFIED AD/»«S PREDI CTOh- CORRECTOR 

C INTEGRATION SCHEME TO SOLVE THE DIFFERENTIAL EQUATION FOR THE 

C STEADY-STATE VELOCITY. 

C 

C P i;s THE VALUE OF THE STEADY- STATE POTENTIAL AT THE STATION# 

C WHERE PREDICTOR-CORRECTOR INTEGRATION COMMENCES; INPUT. 

C DURING THE PROGRAM# F IS CHANGED TO THE VALUE AT CURRENT STATION. 

C H IS THE STEP- size; INPUT THROUGH COMMON BLOCK XI. 

C COMMON BLOCKS XI AND X2 PROVIDE DETAILS OF NOZZLE PROFILE. 

C 

C THE STEADY-STATE QUANTITIES ARE THE OUTPUT# AND 

C ARE PROVIDED BY MEANS OF COMMON BLOCK X5. 

C 

c 

COMMON /XI/ CM# ANGLE# RCC#BCT# GAM# 0#RT#H 
COMMON /XS/ T#R1#R2#NPLAST#NEND#IEXTN 
COMMON /X5/ UC 1000>#DU< IC00)>CClO00;#RWC 1000) 

C 

KP»A 

10 CONTINUE 

PRED e UCNF) + H*<5B.*DU<NP) - S9.*DUCNF-n + 37.*DU<NF-S> 

1 -9.*DUCNP-3>3/2A.0 

P *= P + H 
NP e NP + I 
UP e PRED 

CP *= l.-<GA-M-1.3*UP+.5 

B *5 e>«CF**<-l./ <S.*CGAK-1*D3>3 * C UF**-.25>*A. 

C 

C IF R - 1# THE NOZZLE ENTRANCE HAS BEEN REACHED. 

IF (B-1.3 17# 17# 100 

C 

1? , IF <R-R13 20#2S#2S 

20 DB .= -<<2.*BCT*<R-BT> - < B-RT3*< R-RT) 5> / CRT+RCT-R'3 
60 TO 40 

25 IF (B-B23 30#3S#3S 
30 DR*=-TAN<T3 
GO TO 40 

35 DR » (<2.+RCC*< I .-R3 - < 1 .-E>*< 1 .-R) }**. 53 / tl.-R-RCC) 

40 DO e -.<t)pncDi ,75> « <CP*»!'<<2.^G/^i-13 / < 2.»<GAM- 13) 3 3 / 

1 <0»K1.-CGAM+1.> UP * .53 3 

DUP « DB*DO 

COR = U(NP-13+H* <9 .*DUF+ 19.*EU<NF- 13 - 5.*DU(NF-£) 

2 +DVCKP-3) 3/24.0 

UP = C251.+C0R + 19.+PRED) / 270. 

CP t l.-tGAM-l.l+UP+.S 

B * Q*<CP»*<-1./ <2.*<GAM-I.3>33 * <UP**-.25>*4« 

C 


c 

IF 

R 1# 

THE NOZZLE ENTRANCE HAS BEEN REACHED 


IF 

< R- 1 .3 

62# 62# 100 

U 

62 

IF 

(R-Rl) 

65#70#70 

65 

DR 

= -CC2 

.♦BCT*<R-BT3 - <B-BT>*(R-RT3 >**.S> / <BT+ECT-R3 


GO 

TO 8 5 


70 

IF 

(R-R8) 

7 5#80#60 


ORIGINAX, PAGS 3^ 

OF POOR QUALITS! 

59 - 



75 - DR c *5TPN<T3 

GO TO 85 

60 DR = <f 2.*RCC*« 1 --R) ^ ( 1 .-R>*< 1*-R) >**.5) / (l.-R-RCC)' 

85 DO * -CUF+’t'.TB) * CCP**t <2*+GfiK- 1) / CS.^CGiy*!- 1) > ) ) / 

‘ 1 CQ*< I .“(GAM+ U) * UP ♦ «5)) 

IF CNP .GT« 1000> GO TO 87 
C ' 

C STORE STEADY STATE ODANTITIES AT STATION NP IN RESPECTIVE ARRAYS. 

DU<NP)t=DR*DO 
U(NP> B UP 
C<NP> CP 
EW<KP> »= R 
C 

87 GO TO 10 
100 KPLAST= NP- 1 
RETURN 
END 
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SUBROUTINE COEFFS 03# DU# CO 
C 

C THIS SUBROUTINE COMPUTES THE COEFFI Cl ENTS* 

C U#DU#C#R ARE THE STEADT-STATE OUANTITIES AT THE AXIAL LOCATION# 

C WHERE THE COEFFICIENTS ARE REQUIRED. 

C CC ABE THE COMPLEX COEFFICIENTS. 

C SUBROUTINE INT6BL PROVIDES ALPHA & BETA# THE VALUES OF TRANSVERSE 

C integrals THROUGH COMMON BLOCK X7. 

C 

COMMON /X3/ VC# SVN# 1 P#MODE# NU»KPf 3) 

C0MM0N/X7/ ALPHA! 5# 35# BETA! 9> 35 
COMPLEX CC(255 
DATA GAM /"I. 2/ 

C 

GMINl c GAM - 1. 

M = MODE 

AAB6 «= ALPHA <4#M> * BETA C6#M) 

RSQR = R * R 
C 

LINEAR COEFFICIENTS ****♦♦«****♦>(<#♦*♦♦♦♦♦*’»=*♦♦♦*♦♦♦*♦*♦ 

C 

CCB «= U * <C-U> 

CCCn t CMPLXt CCB# 0-05 
CCR e - U+DU / C 
CCC25 = CMPLXC CCR# 0.05 

CCR *: C ♦ CBETA C8#M5 - BETA <7#M>> / CRSQR * BETA <6#M5 5 
CC<35 »= CMPLXt CCR# 0.05 

CCR c 2. ♦ C ♦ BETA <7#M) / CRSQR * BETA C6#M55 
CCC45 e CMPLXt CCR# 0*0 5 
C 

CCR c C * ALPHA t5#M5 ♦ BETA C9#M5 / tRSQE ♦ A4B65 
CC<55 «= CMPLXt CCR# 0« 03 
CCR « 0*0 

CCI e - 2. * VC ♦ U * KP<M5 
CC16> e CMPLX tCCR#CCI5 
CCR ^ 0.0 

CCI « - GMINl * VC * KPtM) ♦ U * DU / t2« * C5 
CCC75 c CMPLX t CCR# CCI) 

CCR = tVC * KPtM >5 
CCI * 0.0 

CCC8) c CMPLX t CCB# CCI) 

IF tip .NE* 35 GO TO 110 

Q 

C********** NONLINEAR COEFFICIENTS **<^******************************* 
C 

A1 = ALPHA tl»M) 

A2 ALPHA t2#M5 
A3 = ALPHA t3»M) 

B1 *= BETA tl#M) 

B2 >= BETA C2#M5 
B3 « BETA C3#M> 

B4 = BETA C4#K5 
B5 c BETA C5#M> 

CCR - .5 Al+Bl * VC+U / A4B6 
CCI fc CCR 

CCt9) = CMPLX t CCR# CCI 5 
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c 


c 


c 


no 


CCR = - .5 * AI ♦ B3 ♦ VC / tBSQB ♦ A4B6) 

CCI = CCR 

CCUO) e CMPLX (CCR# CCI > 

CCR = - ,5 ♦ AS*BS * VC ✓ (RSOR ’l' A4B63 
CCI c CCR 

CCdl) e CMPLX < CCR# CCI) 

CCR = - ((GAM+IO * U*U ♦ / C4*+3» 1A159S7+A4B6) 

CCI *= - CCR 

CCdS) = CMFLX < CCR# CCI) 

CCR = - (U * A1 ♦ B3) / <4. ♦ RSOB ♦ Aj!iB6) 

CCI »= - CCB 

CC<13> <= CMPLX (CCR# CCI) 

CCR *= - <U ♦ AS ♦ BS) / <4. ♦ RSCR ♦ AAB6i 
CCI * - CCR 

CC(14) = CKPLX <CCR#CCI) ^ 

CCR e “ 3*+U « <1. + .5*GMIM * LM=DU/C) + Al+Bl / <8.*AAB6) 
CCI c - CCR 

CC<15) t CMPLX <CCR#CCI) 


CCR c - DU <1* - (2. -GAM) * U/C) * A1 ♦ B3 / (16 * RSQR * A4B6) 
CCI - CCR ' 

CCC16) = CMFLX (CCR# CCI) 

CCB = - DU ♦ (1. - (2*-GAM) * U/C) * AS * B2 / (16 * RSCR * AAB6) 
CCI = - CCR 

CC(17) = CMPLX (CCR# CCI) 

CCR = - (GKINl ♦ VC + A1 * El) / (A» * AAE6> 

CCI = CCR 

CC(18) “ CMFLX (CCR# CCI) 

CCR e - (GMINl ♦ VC + U + DU ♦ A1 * El) / (A. * C * A-ftB6) 

CCI *= CCR 

CC(19) = CMPLX (CCR# CCI) 

CCR = - GKINl * VC * Al' ♦ - E5> / <4« * RSQR'* AAE6) 

CCI = CCR 

CC(80) = CMFLX (CCR# CCI) 


CCR = - GMJNl * Al * B5 / (2. * RSCR ♦ A/>E6) 

CCI = CCR 

CCCSl) « CMPLX (CCR# CCI) 

CCR s - GMINI ♦ A3 * B2 / (A. * RSQR ♦ A4E6) 

CCI = CCR 

CC(S2) *= CMPLX (CCR# CCD 

CCR s - GKINl ♦ U*A1 * (BA - B5) / (4. * RSCR * AAE6) 
CCI c - CCR 

CC(23) = CMPLX (CCR# CCI) 

CCR = - GKINl + U + Al + B5 / (2»*RSQR * AAB6> 

CCI = - CCR 

CCC8A) = CMPLX (CCE#CCI> 

CCR c - GMIMl * U * A3 * B2 / (A.*R£OR * AAE6) 

CCI = - CCFs 

CCCSS) = CMFLX (CCR# CCD 

CONTINUE 

RETURN 

END 
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SUBROUTINE IKTGRL 

THIS SUBROUTINE EV/iLUATES THE DIFFERENT TRAN.SVERSE INTEGHALS* 

C0KM0N/X7/ ALFHAtS^Si^ BETA(9>3) 

51 1.84118 

52 ® 3.05484 

53 *: 3.83171 
PI * 3*14159 27 

♦*************TANGENTIAL INTEGRALS+*»****>»**+***-»'*-*-************** 

DO 20 NOPT *= 1*3 
£0 ALPHA CN0PT^1> «=0* 

ALPHA C4.1) *= 1*0 

ALPHA *= -1*0 

ALPHA <1*S) = 0*5 

ALPHA t8#2> c -0.5 
ALPHA <3>2> *= -0.5 

ALPHA <4* 2) >= 1.0 

ALPHA <5*25 ** -4*0 

ALPHA <1*35 *= 1*0 

ALPHA <2*35 *= 1*0 

ALPHA <3*35 = -1*0 

ALPHA <4*35 « 2*0 

ALPHA <5*35 » 0*0 

DO 30 -1 = 1*5 
DO 30 U = 1*3 

30 ALFKA<1*U5 = PI+ALPHA<I*U5 

C ♦♦♦♦♦♦♦♦♦♦♦♦♦RADIAL inTEGRALS^4^^4^******************’*‘************ 

C 

DO 40 MODE «= 1*3 
GO TO ( 110* 120* 1305* MODE 
110 M=l 
S*=S1 

GO TO 140 
120 M *=2 

s=sa 

GO TO 140 
130 M*=0 

S=S3 

140 CONTINUE c» 

BETA <1*M0DE> *= RAD2 < 1* 1* 1*M* SI* SI* 5* 

BETA <2*M0DE5 = RAD2 <2* I* I*M* SI* SI* S> 

BETA <3*MODE5 = BAD2 < 7* l> 1*M* SI* SI* S5 
BETA <4*MODE5 = RAD2 <8* 1* 1*M* SI* SI* S5 
BETA <5*MODE5 = RAD2 < 5* 1* 1#M* SI* SI* S5 
BETA (6*«0DE5 = RADI <1*M*S5 

beta <7*M0DE) “ BADl <4*M*S5 

BETA <8*M0DZ) *= RADI <5*M*S5 

BETA (9/M0DE5 «= RADI (8*M*S5 

AO CONTINUE 
RETURN 
END 
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FUNCTION HADl (KOPT^M^B) 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 


c 


c 

c 

c 


c 

c 


c 

c 


115 


1 17 


THIS SUBROUTINE CALCUL/ITES THE INTEGRAL OVER THE INTERVAL 
<0#1) OF THE FOLLOWING PRODUCTS OF TWO BESSEL FUNCTIONS 


KOPT *= 1 UMCB%R> ♦ UM<E*R) + R 
NOPT = S JK(&5:R) * JK(B+R)/R 


NOPT e 3 UFK(B*R) * UMfB*R) * R 
NOPT = A OFM<B+R) * 

NOPT c 5 JPPM(B*H> ♦ JM(B+R> * H 

UM IS THE BESSEL FUNCTION OF FIRST KIND OF ORDER M 
JPM IS THE DERIVATIVE OF JN WITH RESPECT TO R 
JFPM IS THE SECOND DERIVATIVE OF OM WITH RESPECT TO R 
N IS A NON-NEGATIVE INTEGER 
B IS A REAL NtWBER 


DIMENSION FUNCT«SOO> 

DOUBLE PRECISION m* DH> DSTEP^ DR, ARG> EESl^ BES2# EESH, 
1 BESL> PROD> FUNCT* £1/ SS# S3 


NN * 100 
DN = NN 
DH = 1*0/DN 
NPl K NW + 1 

CALCULATION OF INTEGRANDS 

DO 160 1 = 1* NPl 
DSTEP =1-1 
DR = DH * DSTEP 
AEG = B + DR 

CALCULATE BESSEL FUNCTIONS. 

CALL UBES(M*ARG*BES2> S500) 

BESl = BESS 

IF (NOPT »LT. 33 GO TO 130 

CALCULATE FIRST DERIVATIVES OF BESSEL FUNCTIONS* 
CALL JBE£(K+1,ARG>BESH*$5003 
IF (NOPT *E0* 5) GO TO ISO 
IF (I .Efi. I> GO TO 115 
RM = M 

BESl = B * <EM*BES1/AKG - BESH3 
60 TO 130 

IF CM *E0. 03 GO TO 117 
CALL JBESCM- 1* ARG*BESL* S5003 
BESl = B * (BESL - B£SH3/S*0 
GO TO 130 

CALL JBES( 1> ARG*BES1* S5003 
BESl - -BESl * B 
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eo TO 130 

CALCULATE SECOND DERIVATIVES OF BESSEL FUNCTIONS* 

120 IF Cl *£0. 1) GO TO 122 

F « RM * crai * l*0>/’<ARe * ABG> 

EESl = CCF • l.Ol * BESl + EESH/ARG) * B * B 
GO TO 130 

122 CALL JBESCM+2* ABG^BESH^ $500) 

IF CM .EG. 0) EESl *= 0*5 * B ♦ B + CBESH - EESU 
IF CM -EQ. 1) BESl = 0»85 * B ♦ B +CBESH - 3.0+BESD 
IF CM .LT. 8) GO TO 130 
CALL JBESCM*.8>ARG*BESL»$500) 

BESl <= 0.25 + B ♦ B ♦ CBESL - 2.0+BESl + BESH> 

130 FEOD B BESl ♦ BES2 

CALCULATE WEIGHTING FUNCTIONS AND LIMITS FOR R «= 0* 

IF CNOFT .EQ. 2) GO TO lAO 
IF CNOFT *EQ. A) GO TO 150 
FUNCTCI) = PROD + DR 
60 TO 160 

140 IF Cl .EQ. 1) GO TO IA5 
FUNCTCI) = PBOD/DR 
GO TO 160 
IAS FUNCtU) *= 0*0 
GO TO 160 

150 FUNCTCI) = PROD 
160 CONTINUE 


t#*t******^*** SIMPSONS RULE INTEGRATION <^**^************^******** 
NMl = NN - I 

51 c FUNCTCI) + FUNCTCNPl) 

58 e 0.0 

S3 = 0.0 

DO 20 I = 2j NN> 2 

52 » S2 + FUNCTCI) 

20 CONTINUE 

DO 30 I *= 3j NMl# 2 

53 = S3 + FUNCTCI) 

30 CONTINUE 

RESULT = IB ♦ CSl + 4.0>^S2 + 2.0*S3)/3«0 
RADI <= RESULT 
60 TO 501 

500 WRITE C6# 6000) 

6000 FORMAT ( IHl# lOHERROR JBES) 

501 CONTINUE 
RETURN 
END 
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FUNCTION RADS CNOPT#L^W,N,i^^ B/O 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


THIS SUBROUTINE CALCULATES THE INTEGRAL OVER THE INTERVAL 
CO>n OF THE FOLLOWING PRODUCTS OF THREE BESSEL FUNCTIONS 


NOPT = 1 
NOPT 2 
NOPT = 3 
NOPT *= 
NOPT «= 5 
NOPT = 6 


0LCA*R3 * iJW<B+E5 ♦ JN(C+R3 ♦ R 
JL<A*R) * JM<E*E3 * UN<C*R3/R 
OL(A*R3 * UMCB+R) * JN< C*^R3 /CR*R) 
JPLCA+R3 >* ON<B’S'R> « JNCC-t'RJ >♦' R 
JPL(A*R3 ♦ JM(B+E3 * JN(C*H3 
JFL<A+R) * JMCB^R) ♦ dN<C*R>/R 


NOPT B 7 JPL(A*R> * dPMCB+R) + UN<C't‘R).* R 

NOPT B 8 dPPL(A*R5 * JMCB*R) * JN(C*R) + R 

NOPT = 9 JPP1.(A*R3 * dPMtB+R3 * JNCC+RJ * R 

UL IS THE BESSEL FUNCTION OF FIRST KIND OF ORDER L 
dPL IS THE DERIVATIVE OF JL WITH RESPECT TO fi 
JPPL IS THE SECOND DERIVATIVE OF dL WITH RESPECT TO B 
L^ N# N ARE NON-NEGATIVE INTEGERS 
A> C ARE REAL NUMBERS 


DIMENSION FUN CT< 2003 

DOUBLE PRECISION DN# DH* PSTEP# DR> ARGl^ ARG2> AhG3> 

1 BESl^ BES2> BES3» BESH^ BESLj PROD* 

2 FUNCT> BESL1M> SI, S2, S3 

NN = 100 
I»J B NN 
DH = I.O/DN 
NPl B NN + 1 

CALCULATION OF INTEGRANDS **»♦♦♦>!=***♦♦♦♦♦♦♦*=(=«#***+♦ 

DO 160 I B 1, NPl 
DSTEP B I - 1 
DR B DH DSTEP 
AB61 = A « DR 

ARG2 B B * DR 

ARG3 B c DR 


CALCULATE BESSEL FUNCTIONS. 

CALL JBES(N, ARG3jBES 3> S5003 
CALL *JBESCL,ARG1,BES1, £5003 
CALL dBESCM>AF(G2, BESS, £5003 

IF tCNOPT *E0. 73 .OR. {NOPT . EO. 933 GO TO 105 
GO TO 110 
C 
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C CALCULATE FIRST DERIVATIVES OF BESSEL FUNCTIONS. 

105 CALL JBES<M+1jARG 2>EESH* S5003 
IF Cl «E0. 1) GO TO 107 
BM «= M 

BFS2 B B ♦ CFW*BES&/ARG£ - BESH3 
GO TO 110 

107 IF CM .EO. 0) GO TO 109 

CALL JEESCM- 1>ARG2>BESL> S500) 

BF.se «= B + CBESL - BESH)/2.0 

GO TO no 

109 CALL UBESC nA.RG2#BESa> T500> 

BES2 *= -BESS ♦ B 

110 IF CNOFT .LT. A) GO TO 130 
CALL JPESCL+ n ARGnSESHi S500> 

IF CNOPT .GT. 7) GO TO 120 
IF Cl .EO. 1) GO TO 115 

BL *= L 

BESl s= A * CRL*BES1/ARG1 - EESH> 

GO TO 130 

115 IF CL .EC. O) ‘ CO TO 117 

CALL JBESCL- 1, APG1»BE£L, S500> 

BESl = A * CBESL - BESHJ/2.0 
GO TO 130 

117 CALL JBESC l.AROnBESn S500J 
BESl = -BESl ♦ A 
GO TO 130 
C 

C CALCULATE- SECOND DERIVATIVES OF BESSEL FUNCTIONS. 

120 IF Cl .EC. 1) GO TO 122 
RL = L 

F = RL * CRL - NOJ/CARGI + ARGl) 

BESl «= CCF - 1.0) * BESl + BESH/ARGl) * A + A 
GO TO 130 

122 CALL JEESCL+2# ARGn BESH# 5500) 

.IF CL .EO. 0) BESl = 0.5*A*A* (EESH - BESl> 

IF CL .EO* 1) BESl e 0.25 * A * A Jf=CBESH - 3«0*BES1) 
IF CL .LT. 2) 60 TO ISO 

CALL JBESCL-2>ARG1>BESL^S500) 

BESl = 0.25 * A A + CBESL - 2.0=t=BESl + EESH) 

C 

130 PROD c BESl » BES2 * BESS 
C 

C CALCULATE WEIGHTING FUNCTIONS AND LIMITS FOR B «= 0- 
IF CCNOFT .EG. 2) .OR. CNOPT *EQ. 6)> GO TO 133 
IF CNOPT .EQ. 3) GO TO 136 
IF CNOPT .EO. 5) GO TO lAO 
FUNCTCI) tr PROD * DR 
GO TO 160 

133 IF Cl .EQ. 1) GO TO 13A 
FUNCTCI) c PROD/DR 

GO TO 160 

134 BESLIM B 0.0 

IF CNOPT oEO* 6) GO TO 135 


IF 

CCL.EQ. I) 

.AND. 

CM.EG.O) 

.AND. 

CN.EQ.O)) 

BESLIM « 

A/2.0 

IF 

CCL.EQ.O) 

.AND. 

CM.EO. 1). 

.AND. 

CN.EQ.O)) 

BESLIM B 

B/2.0 

IF 

CCL.EG.O) 

.AND. 

CM.EQ.O) 

• AND* 

CN.EQ. 1>> 

BESLIM » 

C/2.0 
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GO It) 155 

135 IF C<L.EQ«0) .AND. CM«EQ»0) 
IF C<L.E0.1> .AND. <MiEG. U 
IF t<L.E6.1> .AND. CM.EO.O 
IF <<L.EQ.£> .AND. CM.EQ.O) 
GO TO IBS 

136 IF (I .EG. n eO TO 13B 
FUNCTtl) « PB0D/CDB*DR) 

GO TO 160 


136 

BESLIM = O.G 



IF 

<<L.EQ.2> 

-AND. 


IF 

tCL.EG.O) 

.AND. 


IF 

ttL.EQ.O) 

• AND. 


IF 

C(L.EQ.l) 

.AND. 


IF 

I (D.EQ. 1> 

.AND. 


IF 

<<L.EQ.O> 

.AND. 


GO 

TO 155 


140 

FUNCTtl) * PROD 


GO 

TO 160 



1S5 FUNCTtI) » BESLIK 


CN.E0*0> 
<M.EQ.2) 
<M.EQ.05 
(M.EO. n 
CM.EO.O) 
<M.EQ. U 


C 

160 CONTINUE 


-AND. tN.EG.OJ) 
-AND. (N.EQ.O)) 
.AND. <N.E£. 1>) 
•AND* tN.£G.O>> 


.AND* tN.EQ«0>) 
•AND. tN.EO.OI) 
.and. <N.ES.S)> 
.AND* (N.EQ.OI) 
.AND* CN.EQ.UJ 
•AND* (N.EO.l)I 


BESLIM » -MA/S.O 
BESLIM » A+B/A.O 
EESLIK = a^C/A.O 
BESDIM « a^A/A.O 


BESLIM *» A*A/8.0 
BESLIM « B»B/g«0 
BESLIM = C+C/8.0 
BESLIM = A*B/A.O 
BESLIM A+C/A.O 
BESLIM «= B*C/4.0 


♦*:^*«***^****4= SIMPSONS INTEGRATION ♦#*^*=^^**#*****.#^*^***+*’^’^ 


NMl »= NM - 1 

SI FUNCTUT + FUNCT(NFl) 

sa = 0.0 

S3 s= 0.0 

DO 20 I *= S* NNj. 2 
52=52+ FUNCTtll 
20 CONTINUE 

DO 30 I = 3. NMl. 2 
S3 = S3 + FUNCTtl) 

30 CONTINUE . 

result = DH * <Si + 4.0«S2 ♦ 
RAD2 = RESULT 
GO TO 501 

500 WRITE 16 * 6000) 

6000 FORMAT C IHl . 1 OHERROR UBES) 

501 CONTINUE 
RETURN 


2.0*S3>/3.0 
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SUBROUTINE RKTZ Tl> G* DUN, I liK) 


C 

C THIS SUBROUTINE PERFORMS A FOURTH ORDER Kia^GE-KUTTA INTEGHATION 

C TO OBTAIN THE INITIAL VALUES FOR THE PR EDI C TOR- CORRECTOR METHOD* 

C 

C NU IS THE NUMBER OF DIFFERENTIAL EQUATIONS TO EE SOLVED* 

C IF IP = 1j INTEGRATION IS CARRIED OUT FOR ZETA ONLY <NU =2). 

C IF IP *= INTEGRATION IS CARRIED OUT FOR ZETA AND AH <NU = 4)* 

C IF IP = 3^ INTEGRATION IS CARRIED OUT FOR ZETA AND GAMMA (NU =4). 

C IP IS PASSED TO THIS SUBROUTINE THROUGH BLOCK COMMON X3* 

C 

C HIS THE step-size; INPUT* 

C T1 IS THE CURRENT VALUE OF STEADY STATE POTENTIAL; INPUT. 

C G ARE THE VALUES OF THE FUNCTIONS AT THE NEXT STEF; OUTPUT* 

C DUM ARE THE VALUES OF THE DERIVATIVES OF THE FUNCTIONS 

C AT THE NEXT STEP; OUTPUT* 

C DUM ARE OBTAINED BY CALLING SUBROUTINE RKDIF* 

C 

C 

COMMON /X3/ VC#SVN#IP/M0DE>NU^KP(3J 
DIMENSION A(4)>G(4J>GZ<4)^ DUM(4>>FZ(4^ 4) 

AC 1)=0. 

A<2)=*5 
AC 3) = *5 
AC4>=I* 

TZ=T1 

NU=4 

IF CIP.EQ. n KU-2 
DO 10 J=1#NU 
10 GZCJ3=GCJ) 

IK-1 

CALL RKDIFCTZ*GZ>DU-!*IK>IRK> 

DO 25 J=l>NU 
25 FZC l^dJaDUMCd) 

DO 30 1K=2>4 
TZ=T1+AC IK)*H 
DO 35 d=l/NU 

35 GZCd)=GCd3 + ACIK)*H*FZCIK-Ud> 

CALL KKDIFCTZ/GZj DIK'^ IK,IFiK> 

DO 50 d=l»NU 
50 FZC IK,d)=DUMCd> 

30 CONTINUE 

DO 55 d=UNU 

55 GC J)=GC JJ+H’t'CFZC 1, J) + 2*>t-CFZC 2^ J) + FZC 3,d> ) + FZC 4^d) ) /6* 

IK=4 

CALL RKDIFC TZ/G» DUM^ IK/ 1RK> 

75 RETURN 

END 
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SUBROUTINE RKD1F(F,U*GE>IK« 1EK> 


THIS SUBROUTINE EVALUATES THE DIFFERENTIAL ELEMENT IN THE 
RUNGE-KUTTA INTEGRATION SCHEME. 

FIS THE CURRENT VALUE OF STEADY-STATE POTENTIAL; INPUT- 
G ARE THE VALUES OF THE FUNCTIONS AT P; INPUT- 
OF ARE THE DERIVATIVES OF FUNCTIONS AT PI OUTPUT- 


COMMON /XI/ CM^ANGLEj RCC>RCT>GAK#0#RT^ DP 
COMMON /X8/ T#Ri>R2-NPLAST>NEND- I EXTN 
COMMON /X3/ VC^SVN>IF-K0DE*NU-KP<3) 

COMMON /XA/ RU<7I ^KDUCTJ^ETHRl^GTORl 
COMMON /X6/ AFN-AFNW AFN2 

COMiPLEX AFNC 1000)# AFNK 1000)-AFN2t 1000> 

COMPLEX CCI25)#CFH#CFM#CFN#INHKG 

COMPLEX ZETA>ZETAJ#AH-AH1# CGAM- CG/«^ 1#ZTHHI» 6THR1# AP# API# AF2 
DIMENSION G(4)#GPC4) 

C 

ZR « e< 1) 

ZI = GC2) 

ZETA != CMPLX CZE#Z1) 

GO TO < 110> 120> 130)#IF 
120 AKR = GC3) 

AHI = G<4) 

AH c CMPLX <AHR#AHI) 

GO TO 110 
130 CONTINUE 

GAMR = G<3) 

GAMI B GCA) 

CGAM B CMPLX CCAMR#GAMI) 

110 CONTINUE 

IF CF5 15# 10- 15 
10 GPCl) « REALCZTHRU 
GPC25 B AIMAG(ZTHRl) 

GO TO < lAO# 150# 160)- IP 
ISO AHI = AH * ZETA 

6PC3) = REAL CAHl) 

GPC4) « AIMAG CAHl) 

GO TO 140 
160 CONTINUE 

6PC3) = REAL CGTHRl) 

GPC4) s= AIMAG CGTHRl) 

140 CONTINUE 
GO TO 20 

15 ICL‘=2#ihk-2 

IF CIK -EQ- 1) ICL=2*IRK-3 
IF CJK -EG* 4) ICL B 2+1 RK - 1 
UbRUCI CL) 

DU=BUU< I CL) 

Ct=l,*CCAM-l.)*U*-5 

B=QiM C C)**C - 1/C2*CGAM“ 1 . ) ) ) )*CU*+- * 25)* 4- 
CALL COEFFS C U> DU# C# R# CC) 

CFH = CCC 1 ) 

CFM B CCC2) + CCC6> 

CFN B CCC3> + CCC4) + CCC 5) + CCC7> + CCC8) 
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ZETAl a C -CFM ♦ 2ETA CFN> / CFH * ZETA * ZEtA 
6P<1> *» REAL CZETAl) 

GP<e> a AIMAQ <ZETA1> 

CO TO <170#180»190)^ IP 
180 AHl a AH ♦ ZETA 

GPC3> •= REAL <AH1> 

GP(A> *= AIKAG CAHI) 

60 TO 170 
190 CONTINUE 

GO TO C30* AO>AO» 50>» JK 
30 AP t AFN <lifl<«l> 

API a AFNl <IEK-1> 

APS = AFNS <IBK-1J 
GO TO 60 

AO AP t .5 * (AFN (IRK-1> + AFN <IRK1> 

API ,5 * <AFTM1 tIRK-1) + AFNl <IKKJ> 

APS « ,5 CAFNS CIRK-1> + AFNS <IRK>J 
GO TO 60 

SO AP « Am CIRKJ 

API a AFNl <1RK) 

APS e AFNS tlEKl 
60 CONTINUE 

INHMG « - CC<18) ♦ AP * APS - CC<18) ♦ API * APS - CCCC9> 

1 + CC<15)> ♦ API ♦ API - <CCC13) + CC<14> + CCCI9> 

2 + CC(23) + CCtSAJ + CCC25>) * API * AP - <CC<101 + CCCllJ 

3 ♦ CCC 17J + CCC20) + CCCSn + CCCSS>> ♦ AP ♦ AP 

C6AM1 a € “ ZETA + .5* CGAM-l.) * OU/C - CFM/CFHJ ♦ CGAM 

1 * IKKMG / CC * CFH> 

GPC3J « REAL CCGAMl) 

GPCA) *s AIMAG (CCAMl) 

170 CONTINUE 
20 RETURN 

mu 
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SUBROUTINE ZADAWS (H>X/ Y> BY* I T0K2 ) 

THIS SUBROUTINE CABBIES OUT A MODIFIED ADAMS FREEI CTOE-CORRtCTOR 
INTEGRATION SCHEME TO SOLVE THE VARIOUS DIFFERENTIAL EQUATIONS AS 
DESCRIBED BELOV 

IF IP = 1> INTEGRATION IS CARRIED OUT FOR ZITA ONLYi 

IF IP = Z* INTEGRATION IS CARRIED OUT FOR Z ETA AND AH; 

IF IP *= 3^ INTEGRATION IS CARRIED OUT FOR ZLTA AND GAMMA. 

IP IS PASSED TO THE SUBROUTINE THROUGH COMMON BLOCK X3. 

H IS THE step-size; INPUT- 

X IS THE VALUE OF STEADY-STATE POTENTIAL AT THE STATION ^ 

VHERE THE PREDICTOR-CORRECTOR INTEGRATION STARTS? INPUT* 

DURING THE PROGRAM/ X IS CHANG ID TO VALUE AT CURRENT STATION- 
Y ARE THE VALUES 'AT K / OF THE FUNCTIONS/ VHOSE EQUATIONS ARE 
BEING solved; INPUT AND OUTPUT. 

BY ARE THE BERIVATIUES OF Y; INPUT ANB OUTPUT- 

ITORZ PASSES TO MAIN PROGRAM THE I NFORMATI ON - AS TO VHICH VARIABLE 
(TAU OR ZETA) HAS BEEN INTEGRATED. 

ITORZ = I t INTEGRATION OF EGUATION FOR TAU- 

ITORZ = 2 ; integration of EQUATION FOR ZETA. 


COMMON /XI/ CM/ ANGLE/ RCCjRCT/GAM/Q/RT 
COMMON /XS/ T/KJ/RS/NFLAST/NEND/IEXTN 
COMMON /X3/ VC/ SVN/IF/MO,BE/NU/KP< 3> 

COMMON /X5/ U< 1000>/ DU< 1 0005/ C( 1000)/ RVC 1000) 

COMMON /X6/ AFN/ AFNl^AFNS 
COMMON /X8/ ZETA/ TAU/ CCEXT 
COMPLEX ZETAC 1 000)/ TAUCl 0005/ CCEX TC S5) 

COMPLEX AFN< 1000)/ AFNK 1000)/AFN£C 1000) 

COM PL EX CC<e5)/CFH/CFM/CFN/lNHMG/ZETAl/AH/AHl/AHJ>/Ap/AFl/APS/ 
1 CGAM/ CGAM 1 

DIMENSION Y( 4)/ DYC4/ 4) / EFC 4)/PREDC 4)/ CORI 4) 

NP=4 

ITORZ « 8 

IF <IEXTN -NE. l> GO TO 10 

BEFINE STEADY STATE QUANTITIES IN THE EXTENSION REGION- 

UEXT = UCWEND) 

CEXT - CCNENB) 

BEXT = RWCNEND) 

BUFXT = DUCNENB) 

CALL CO EFFS < UEXT/ BUFXT/ CEXT/ REXT/ CCEXT) 

NU IS THE NUMBER OF EQUATIONS TO BE SOLVED. 

10 CONTINUE 

DO 15 d=I/NU 

FRED( J)=YC J)+H*( 55- + DY< J/ 4) - 59 .#BY C J/ 3)+ 37* =('DY (0/ 2) 

1 -9.*DYCd/ 1 ))/24. 

15 CONTINUE 

X=X+H 
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120 


130 

110 

25 

20 


30 


150 


160 

C 

C 

C 


NF»nP+ 1 
ZR=FBEO< 13 
ZI«=PREDC2> 

EETACNP-3 w CKPLZ <ZB#Z13 
GO TO C110*120»130>j IP 
AHR e PRED<33 
/mi * PRED<4> 

AH « CKPLX </mR>AH1.3 


GO TO 110 
COKTINOE 

CGAK « CHPLX (PHEDC 3)»PRED< A) ) 


C0KTIM3E 

IF (NP ‘EE. NFLAST> GO TO 20 
DO 25 I »= 1>25 
CCtI) = CCEXT(I) 

GO TO 30 
CONTI NUE 
UPelKNP) 

DUP*=DU<NP> 

CP«C<NP> 

ReRW(NP> 

CALL COEFFS f UP#DtJF^CF»R*CC3 

CONTINUE 

CFH *= CCC13 

CFN « CC<2> + CC<63 

rtTM fc CCC3D + CCCA)'+ CC< 53 + CC<73 CCC83 

ZETAl *= < - CFW * ZETA(NP> - CFN) / CFH - ZETA<NF) ♦♦.2 


DP<13 *= REAL <Z ETA 13 
DPC23 « AiWAG <ZETA13 
GO TO < 1A0> 150 j 160)^ IP 
AHl »= AH * ZETA(NF3 
DPI 33 = REAL CAH13 
DP<A3 AIMAG (AH 1 3 
CO TO lAO 
CONTINUE 


AP^APl AND AP2 ABE THE VALUES OF THE AMPLITUDE FUNCTION AND 
THEIR DERIVATIVES AT THE CURRENT STATION. 

AP AFN(NP3 

API AFNKNP3 . _ 

AP2 b AFN2(NP3 


140 


45 


INHMG - - CC( 18> » AP • AF2 - CC( 18> * 

CC(15)> » API • AfJ - CCC<13> * CCOA) * 

CCC23) + CC<243 ♦ CCC25D3 ♦ API ♦ AP - 
rrci73 ♦ CCC203 + CC(813 + CC<2835 ♦ AP » AP 
<- ZETACKP3 + -5* (GAM- 1.3 * DUP/CP - CFM/CFH3 * CG^ 

- INHMG / (CP * CFK3 
(CGAMl 3 
(CGAM13 

CONTINUE' 

^R(U3b Y(U3 + H+(DY<0^23-5.*DT{J^33+19.*ETCJ^43 
I 49 .*DP(J3 3/24.0 

Y(J)« (251.*C0R(J)+19.*PRED(U33A27,0. 


CGAMl 

I 

DP<3) 

DF(43 


4 

4 

4 


REAL 

aim'ag 


CC< I 1) 
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EO 55 IcWNU 
CO 55 

55 * DV<I»J+n 

ZH=Ytl) 

ZIcYCSJ 

ZETACNP7 = CMFLX <ZF.>ZI) 

ZETAl = ( - CFK ♦ ZETACNF) - CFN) / CFH - ZE1ACNF5 **2 
DY C1j*4) *: REAL <ZETA1> 

DY C£,4> e AIKAG <ZETAi) 

GO TO « 170# 180# 190># IR 
IBO AH e CKFLX (YC3J#Y<A)) 

AHl c AH * ZETA<NP1 
DYC3#AJ = REAL <AHU 
DY(A#A} *= Alf^AG <AHI) 

IF CMODE«NE«I) GO TO 182 
AH2 = AHl * 2ETA<NF) + AH ♦ ZETAl 
AFNCNP) *= AH 
AFNKNF) = AHI 
AFN2CNF) n AH2 
188 GO TO 170 
190 CONTINUE 

CGAM = CMPLX CY<3)#Y<A3) 

CGAMl = <- ZETACKF) + . 5* CGAK-l.J * LUP/CF - CFM/CFH> ♦ CGAW 
I - - IKHKG / (CF * CFH> - 

DYC3#A5 «= REAL tCGAMll 
DY (4# A) c AIF/AG (CGAN13 
170 CONTINUE 

IF <NF .EO. NEND) GO 10 100 
X 

C DECIDE tHICH EQUATION IS TO BE INTEGRATED : TAU OR ZETA 

C 

IF (CABS (ZETACNP)) .LT. 103 GO TO 10 
JTORZ = I 
C 

C CALCULATE VALUE OF TAU AND ITS DERIVATIVE AT LAST FOUR STATIONS. 
DO 410 I = 1#4 

410 TAU CNF-4+I) c I */Z ETACNP-4+ 1 3 
YCn *= BEAL (TAUCKP3) 

Y(£3 «= AIWAG (TAU(NP3> 

EO 420 I = 1#4 

TSGR B REAL <TAU(NP-4+I3 * TAUCNF'4+13> 

TSei = AIKAG (TAU(NF“4+I3 * TAUCNP“4+I ) > 

ZPR B DY( 1#I 3 
ZPI B DY(8#I> 

DY(l#t) = - TSCR+ZPR + TSGI+ZPI 
DY(2>I> B - TSOR+ZPl - TSOI +Z PR 
420 CONTINUE 
C 

CALL lADAKS < H# NP# X# Y # DY# I 0* I TORZ 3 
GO TO CIO# 1003# IQ 
100 RETURN 
END 
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SUBROUTINE TARjftMS <H*NF*X*Yjr»V*ie*IT0RZ5 

THIS SUBROUTINE C/^RRII-S OUT A MODIFIED ADAMS FEED! CTOR- CORRECTOR 
INTEGRATION SCHEME TO SOLVE THE VARIOUS DIFFERENTIAL EQUATIONS AS 
DESCRIBED BELOV 

IF IP «= U 1NTECF.ATION IS CARRIED OUT FOR TAU ONLY* 

IF IF «= 2* INTEGRATION IS CARRIED OUT FOR TAU AND AH5 

IF IF e 3i INTEGRATION IS CARRIED OUT FOR TAU AND GAMMA. 

IF IS PASSED TO THE SUBROUTINE THROUGH COMMON BLOCK X3* 

H IS THE STEP-SIZE; INPUT. - 

X IS THE VALUE OF STEADY-STATE POTENTIAL AT THE STATION * 

VHERE THE PREDICTOR-CORRECTOR INTEGRATION STARTS! INPUT. 

DURING THE PROGRAMj X IS CHANGED TO THE VALUE AT CURRENT STATiON. 
Y ARE THE VALUES AT X ^ OF THE FUNCTIONS# WHOSE EQUATIONS ARE . 
BEING solved; INPUT AND OUTPUT. 

DY ARE THE DERIVATIVES OF Y; INPUT AND OUTPUT. 

IQ INDICATES WHETHER INTEGRATION IS COMPLETE; OUTPUT. 

IQ » I 1 INTEGRATION IS TO BE CONIINUEP BY SUBROUTINE ZADAMS. 

IG « 2 * INTEGRATION IS COMPLETE. 

ITORZ INDICATES WHICH EQUATION SHOULD EE INTEGRATED : 

ITORZ 1 J INTEGRATION OF EQUATION FOR ZETA. 

ITORZ » 2 t INTEGRATION OF EQUATION FOR TAU* 


C 

10 

C 


15 


120 


130 


COMMON /XI/ CM# ANGLE# HCC*RCT* GAM# 0#RT 
COMMON /X2/ T#R1#R2#NPLAST#NEND#IEXTN 
COMMON /X3/ VC#SVN#1P#K0DE#NU#KFC3> 

COMMON /X5/ UC1000)#DU< 1000)#C< 1000)#RW( 1000) 

COMMON /X6/ AFNi AFNl# AFN2 

COMMON /X8/ ZETA# TAU# CCEXT 

COMPLEX AFN( 1000)# AFNlt IOOO)#AFN2C 1000) 

COMPLEX CC< 25) # CFH# CFM# CFN# 1 NHMG# AH# AH 1 # AP# AP 1 , AP2# C6AM# CCAK 1 
COMPLEX ZETACIOOO)#TAUC 1000)#TAU1# CCEXTC25) 

DIMENSION Y< A)# DY( A# A)# DP( A)#PRED(A)# CORC A) 


CONTINUE 

NU IS THE NUMBER OF EQUATIONS TO EE- SOLVED. 

DO 15 J »= 1#NU 

FRED< J)=Y< J)+H*<55.*DY< J# A>-59.*DY(0#3) + 37.*DY< J# 2) 
1 -9.*DY(0# 1))/2A. 

CONTINUE 
X = X+H 
NP e NP + 1 
TR B PRED t 1 ) 

T1 * FRED C2> 

TAU <NP) = CMFLX C TR# TI ) 

ZETA CNP) = 1./ TAUCNP) 

GO TO < 1 10# 120# 130)# IP 
-AHR c PRED(3) 

AHI *= PRED <A) 

AH e CMPLX t AHR# AHI) 

GO TO no 

CONTINUE 

CC^ “ CMPLX (FRED(3)#PPED(A>) 
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110 CONTINUE 

IF (NP .LE. NPLAST) GO TO 80 
C 

C OBTAIN COEFFICIENTS IN THE EXTENSION SECTION. 

rO 25 I « 1*25 
85 CCflJ = CCEXTCI) 

C 

GO TO 30 
20 CONTINUE 

DUP = DU(NF) 

UP = UCNF^ 

CP .= C<NF> 

E = RW <NF) 

CALL COEFFS C UP^ LUFj CP* R> CO 
30 CONTINUE 

CFK = CC< 1 ) 

CFF e CC<8) + CCC65 

CFK s CCC3) + CCCAJ + CC( 5) + CCC7) ♦ CCC8) 

TAUl c 1. + (CFM + CFN * TAUCKP3) * TAUCNPJ ✓ CFH 
rP<l) = REAL. fTAUn 
DP<2> >= AINAG < TAUl 5 
GC TO 150-* 1603^ IF 

150 AHl = AH / TAUCNPJ 

RFC 33 *= REAL CAH13 
DFCA) = AIMAG CAHID 
GO TO lAO 
160 CONTINUE 

C 

C AF, API AND AF2 APE THE UALUFS OF THE AMPLITUDE FUNCTION AND 

C THFIR DERIVATIVES AT THE CURRENT STATION. 

AP AENCNP) 

API = AFNHNP? 

APS = AFN2CNP) 

C ' 

INHM6 *: - CCC18> ♦ AP ♦ APS - CCC12) * AP1 * AF2 - <CC<9) 

1 + CCC15)> * API * API <CC<13> + CCCiA? + CC<19) 

2 + CC<23) + CCC24) + CCCS5) ) * API * AP CCCCIO) + CCC 1 IJ 

3 + CCC 17) + CCC 20) + CCC 21) + CCC 22) ) * AP * AF 

CGAKl c < - ZETACNF) + .5 ♦ CGAM - 1.) * DUF/CF - CFM/CFh) * CGAW 

1 - INHMG / CCP * CFH) 

DPC3) .= REAL C CGAKl) 

DPCA> = AIMAG CCGAMl) 
lAO CONTINUE 

DO AS J=1^NU 

COncJ)c YCd) + h*CDYCO»2)-5.*EyCd>3)+19.*EYCO^ A) 

1 +9.+BFCJ) )/2A.O 

AS YCd>= <25l.*C0RCd)+ 19.*FREDCd))/270- 

DO 55 I=UNU 
DO 55 J=l>3 

55 DYCI^d) = DYCI>d-e-l) 

TR = YCD 
TI = YC2) 

T2 »= TF.+ TR + TI*TI 

TAU CNF) = CMPLX <TR#TI) 

ZETA CNP) *= 1./ TAUCNP) 
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J80 


182 

190 


170 

C 

C 

C- 


C 

C 


420 

C 


100 

105 


TAUl «= It + CCFM + CFN * TALtCNP)) * TAUCNF> / CFH 
0Y Cl>4) = PEAL <TAyi> 

DY C2«/|) c AIM AG (TAU1> 

GO TO < 170/ leOi 190)# IP 
AHR o YO) 1 

AKI *= Y<4) 

AH t= CMFLX <AHR/AHI) 

AHl c AH / TAUCNF) 

DY C3/A> = REAL (AHl) 

DY (4/4) =■ AIMAG CARD 
IF (KOLE tKE. 1) GO TO 182 
AFR(NF) = AH 
AFNKNF) *= AHl 

AFN2 (NF) c ( TAIKNF) * AFWKMF) - TAUl ♦ AFN<fr)f) ) / 
i < TAUCKF)tTAU(KF) ) 

GO TO 17 0 
CONTINUE 

CGAM = CKFLX (Y(3)/Y(4)) 

C6AM1 = ( - 2ETACNP) + •? * <GA« - It) ♦ LUP/CP - CFM/CFH) 
I - INHKG / C CP ♦ CFH) 

DV (3/4) = REAL (C6AN1) 

DY (4/4) *: AIKAG (CGANl) 

CONTINUE 

IF <KP .EQt NEND GO TO 100 

decide VHICH EGUATION is to PE integrated : TAU OR 

IF <CABS(TAU(NP)) *LT. 10) GO TO 10 
ITORZ = g 

Y(l) c REAL < ZETACNPJ ) 
y<2) t: AIMAG ( ZETA<NP) ) 

calculate DERIVATIVES OF ZETA AT THE LAST FOUR POlNTEt 
DO 420 1 *= 1/4 

ZSQR = REAL < ZETAtKp-4+lJ ♦ ZETA(NF-4+I) )’ 

ZSQl «= AIKAG C ,ZETA(NP-4+1 ) * ZETA<NF-4+I’') )' 

TPB •= 

TPI = DYC2/I) 

DYC1>1) = - ZSQR+TPR + ZSGI*TPI 
DYC2/I) *= - ZSGR+TPi - ZS61 + TPR 
CONTINUE 

10 = 1 ! ' 

RETURN 
IQ »= 8 
RETURN 
END 


* CGAK 


ZETA 
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APPENDIX B 

PR0GEAM.C0EPFS3D: A USER'S. MANUA: 

. i . :■.! \ ... 

Program COEFFS3D calculates the coefficients of both .the linear and 

V - 

nonlinear terms that appear in Eq_. (20) . These coefficients are required as 
input for Program LCYC3D (see Appendix c) which numerically integrates this 
system of equations. Program C0EFFS3D is a slightly modified version of the 
program described in detail in Appendix C of Ref. 11. The modification lies in 
the evaluation of one more coefficientj p) defined by 






e 


r 

p 



® ®!d0 
P 3 



R R.rdr. 
P 3 . 


This coefficient represents the effect of nozzle nonlinearities. Except 
for this additional coefficient, the two programs are very similar in the 
structure of their numerical calculations and their output. Hence in this 
user's manual, only the listing of the entire program together with a precise 
description of the necessary input -is given. For details of the program, 
one is referred to Appendix C of Ref. 11. 

In the following description of the input, the location number refers 
to columns of the card. Three formats are used for input: "A" indicates 

alphanumeric characters,- ”l"‘ indicates integers and "F" indicates’ real 
numbers with a decimal point; -For the "l" and '"F" formats the values are 
placed in fields of five and ten locations respectively (right justified). 


Ho. of 


Cards 

Location 

Type 

input Item 

Comments 

1 

1-72 

A 

Title 

Title of the case 

1 

1-10 

F 

GAMMA 

Ratio of specific heats 


■ 11-20 

F 

UE 

Steady-state Mach mmiber 
at nozzle entrance 


21-30 

F 

RED 

Length-to-diameter ratio 


31-ij-O 

F 

ZCOMB 

Length of the combustion 
zone 
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Location 

Tyne 

Inpnt Item 

Comments 

lj.l-.45 

I 

KLROPS 

If 0: droplet momentum 
source neglected 

If 1: droplet momentum 
source included 

46-50 

I 

5F0ZZIE 

If 0: quasi-steady nozzle 
If 1: conventional nozzle 

1-5 

I 

Njmx 

Humber of series terms 
(complex) 

6-10 

I 

KOULIM 

If 0: linear terms only 
If 1: both linear and 
nonlinear terms 

11-15 

I 

NEGL 

If 0: all non-zero coeffi- 
cients calculated 
If 1: small coefficients 
neglected 

16-20 

I 

mm 

If 0 : printed output only 
If 1: printed and ’written 
into file 

If 2: -written into file only 
If 3 : card output orOy 

21-25 

X 

lOZULl 

If 0: nozzle nonlinearities 
neglected 

If 1: nozzle nonlinearities 
included 

26-30 

I 

DTZDATA 

If 0: nozzle admittance values 
input through cards 
If 1: nozzle admittance 
values input through file 
If HZDATA is 1, ITOUI in 
program UOZADM should be 1 

card is necessary only if HEGL - 1. 

* 
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Ho. of 
Cards 

Location 

Type 

Input Item 

Comments 

1 

1-10 

F 

SMI 

Linear coe_fficients with 
absolute value less than' 
SMI neglected 


11-20 

F 

SJVI 2 

nonlinear coefficients 
■with absolute value less 
than SM 2 neglected 

The next' 

HJMAX cards are 

necessary 

only if NOZZIE 

= 1 and HZDATA /=. 0 . 

HJMAX 

1-5 

I 

J 

Integer which, identifies 
the series term 


6-15 

F 

AMPL(J) 

Amplitude of the linear 
nozzle admittance 


16-25 

F 

PHASE (J) 

Phase of the linear nozzle 
admittance 

The next 

HJMAX cards are 

necessary oxQy if KZOATA 

= 0 and NOZJTLl = 1 . 

HJMAX 

1-5 

■I 

J 

Integer which identifies 
the series term 


6-15 

F 

GHOZ(J) 

Real part of the nonlinear 
nozzle admittance 

- 

16-25 

F 

GHOZ(J) 

Imaginary part of the 
nonlinear nozzle admittance 

HJMAX ■ 

1-5 

■ I 

J 

Integer which identifies 
series term 


6-10 

I 

L(J) 

Axial mode number, ^ 


11-15 

I 

M(J) 

Tangential mode nuniber, m 


■ 16-20 - 

I 

N(l) 

Radial mode number, n 


21-25 

I 

HS(l) 

Urs(l) = 1 : ®. = sin (mB) 
HS ( I) = 2 : = cos (me ) 


26-30 

A 

HAME(I) 

Four character name 
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FORTRAN Listing 


c 

c 

C pJtOGpAM COEFFS3D ♦♦*♦*♦*♦+*+*♦*♦++*«♦+♦♦♦♦ 

C 

C THIS PROGRiiK COMPUTES THE COEFFICIENTS WHICH APPEIAh 

C IN THE differential EBUATIONS WHICH GOVERN THE MODE- AKFLI TUDE 

C FUNCTIONS. THESE COEFFICIENTS AHE PWCHED ONTO CARDS FOB 
C INPUT INTO PROGRAM LIMCYCa 

C 

C THE FOLLOWING INPUTS ARE REOUIREDJ 

C THE TITLE OF THE CASE® 

C GAMMA IS THE SPECIFIC HEAT RATIO. 

C UE IS THE STEADY STATE MACH NUMBER AT THE NOZZLE ENTRANCE. 

C RLD IS THE LENGTH-fO-DI AMETER RATIO. 

C ZCCMB IS THE LENGTH OF THE HFOION OF UNIFORMLY DISTRIBUTED 

C COMBUSTION. EXPRESSED AS A FRACTION OF THE CHAMBER LENGTH. 

C . NDROFS DETERMINES THE PhESEFJCE OF DROPLET MOMENTUM SOURCES* 

C NDROPS = 0 DROPLET MOMENTUM SOURCE NEGLECTED. 

C NDROPS « I DROPLET MOMENTUM SOURCE INCLUDED. 

C NOZZLE SPECIFIES THE TYPE OF NOZZLE USED* 

C NOZZLE *= 0 QUASI -STEALV 

C NOZZLE fc 1 CONVENTIONAL NOZZLE. 

C -FOR CONVENTIONAL NOZZLE 

C AMPL IS THE NOZZLE AMPLITUDE RATIO. 

C PHASE IS THE NOZZLE PHASE SHIFT. 

C NOZNLl DETERF/INES THE PRESENCE OF NOZZLE NONLINEARITI EE 
C NOZNLl *= 0 NOZZLE NON-LINEARITIES NEGLECTED. 

C NOZNLl t 1 NOZZLE NONLINEARITIES INCLUDED. 

C NZDATA DETERMINES HOW THE NOZZLE DATA IS SUPPLIED 

C NZDATA = 0 FROM CARDS* 

C NZDATA fc 1 FROM A FAETHAND FILE. 

C NJMAX IS THE NUMBER OF KOLE-AMPLl TUDE FLtOCTIONS IN THE ASSUMED 

C SERIES SOLUTION. NdMAX MUST NOT EXCEED MX. 

C THE COEFFICIENTS COMPUTED ARE DETEFMINED BY NOKLIK AS FOLLOWS 

C NONLIN c 0 LINEAR COEFFICIENTS ONLY 

C NONLIN = 1 BOTH LINEAR AND NONLINEAR COEFFICIENTS 

C COEFFICIENTS TO BE NEGLECTED APE DETERMINED BY NEGL 
C AS FOLLOWS: 

C NEGL »= 0 TERMS SMALLER THAN 0-00001 ARE NEGLECTED. 

C NEGL = 1 LINEAR TERMS SMALLER THAN SMI AND NONLINEAR 

C TERMS aOALLER THAN SMg AHE NEGLECTED* 

C THE OUTPUT IS DETERMINED BY NOUT AS FOLLOWS 

C NOUT « 0 PRINTED OUTPUT ONLY 

C NOUT = 1 PRINTED AND WRITTEN ON FASTRAND FILE. 

C NOUT' = 8 FASTRAND FILE ONLY* 

C NOUT = 3 CARD OUTPUT ONLY- 

C ■ EACH MODE- AMPLITUDE IS ASSIGNED AN INTEGER J. 

C THE MODE IS SPECIFIED BY THE INDICES L<J>. K<d). AND N<J5* 

C L<J> IS THE AXIAL MODE NUMBER AND MUST NOT EXCEED 5. 

C MCJ) IS THE AZIMUTHAL MODE NUMBER AND MUST NOT EXCEED 8* 

C NCJ> IS THE RADIAL MODE NUMBER AND MUST NOT EXCEED 5. 

C THE INTEGER NSCd> IS ASSIGNED AS FOLLOWS* 

C NS e I A-FUNCTION £INCM>nHETA) + COSH<I + B*Z) 

C NS - 8 B- FUNCTION COSCM*THETA> * CO£HCI*E+Z> 

C '• NAMECd) IS A FOUR-CHARACTER NAME. 
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PAHAKETEF 

DlMEKSIOtj 


COMMON 


MX=5^ MXSsIOj MX4m=£0 

LCMXJj NCMXJj NAMEtMX)* SCMXJ> SJ(MX)i TITLEC80^j 
RJEOOT< 10» SJVALUO^S)# Cl<MX2>MX2>j CC *>MX2^ KX2)# 
DCKXa,KX2>MX2>> AMFL<MX)/ PHASECMX)# AZK2), 
EEE1C9^9^93» BES2C9^9j9)> BES3<9>9>9>* 

VC25> JCCMX2)-» TSC/i*MX2)> T£6<KX2)> KMAX<5) 

CRSLT> CI> ZEJ> ZEP1» ZEF2# CZEj CAZ* CRAD« 

Gl> BCOEFj C6AM# CAX^ ECMX>^ BC(MX>> YKOZCKXJ^ 
CNORKCKX)^ C££0<MX)* TAKINTC25j RADINT(3)^ 

AXINTC/|»3>» CC<5#MX>MX># CD! (MX^MX# MX>* 

CD2XKX,MX^MX>> AXf45- Tl> T2# Dl» D2* D3^ 
CD3<MX>MX>MX>, CD4 CMXjMX»HX>j» GN0Z<MX5 
B /BLK2/ M<MX>^ NSCKX> 


BATA INPUT* 

FI - 3.1/115927 
£Wl c 0.00001 
SWa n 0*00001 
SM3 *= O.OOOOl 
Cl <= C0*0*1.0> 



INPUT ROOTS 

AND VALUES 

OF BESSEL 

. FUNCTIONS. 


DATA (CFJEODTCl/ J)/ J 

= 1/5)/ 1 

= 1/9)/ 


1 

3.8317 1/ 

7*01559/ 

10. 17347/ 

13*32369/ 

16.47063/ 

2 

1.641 18/ 

5.33144/ 

8.53632/ 

11*70600/ 

14.8 6 3 59/ 

3 

3.05424/ 

6*70613/ 

9.96947/ 

13* 17037/ 

16.34752/ 

4 

4.201 19/ 

8*01524/ 

1 1.34598/ 

14. 58 585/ 

17*78875/ 

5 

5.31755/ 

9.28240/ 

12.68191/ 

1 5*9641 1/ 

19. 19 603/ 

6 

6.41562/ 

10.51986/ 

13.987 19/ 

17.31884/ 

20* 57551/ 

7 

7.50127/ 

1 1*7 349 4/ 

15.26818/ 

18*63744/ 

21.93172/ 

6 

6.57784/ 

12*93839/ 

16. 52937/ 

19.94165/ 

23* 26805/ 

9 

9.64742/ 

14.11558/ 

17.77401/ 

£1.22906/ 

24.56 7 80/ 

DATA <CHJVAL(I/ J>/ J = 

1/5)/ 1 e 

1/9)/ 


1 

-0.40276/ 

0.30012/ 

-0.24970/ 

0.21836/ 

-0* 19647/ 

2 

0.58187/ 

-0.34613/ 

0.27330/ 

-0.23330/ 

0*20701/ 

3 

0*48650/ 

-0.31353/ 

0*25474/ 

-0.22068/ 

0* 19794/ 

4 

0*43439/ 

-0.29116/ 

0.24074/ 

-0.21097/ 

0* 19 048/ 

5 

0*39965/ 

-0.87438/ 

0.889 59/ 

-0.80276/ 

0* 18 403/ 

6 

0.37409/ 

-0.26109/ 

0.88039/ 

-0. 19 560/ 

0* 17649/ 

7 

0.35414/ 

-0.25017/ 

0.81861/ 

-0. 189 78/ 

0. 17363/ 

8 

0*33793/ 

-0.24096/ 

0.20588/ 

-0. 18 449/ 

0. 169 29/ 

9 

0*32438/ 

-0.83303/ 

0. 19998/ 

-0. 17979/ 

0. 16539/ 


INPUT PAh'ANETEBS. 

BEAD (5>5000> END = 600> <TITLE<I>a 1 li* 72) 

BEAD (5*. £001) GANMA> UE, KLD> ZCOMB* NtROPE^ NOZZLE 
IF-CGAKKA) 600/ 600/ 8 

HEAD C 5/ 5004) NWjAXj N0NLIK> NEGL* NOUT/ NOZNLl/ NZDATA 
IF CNEGL .EG* 1) READ ( 5/5005) St-ll/ £N2 
IF (NOZZLE *E&» l> GO TO 5 
COMPUTE AENITTANCE FOR QUASI -STFADY NOZZLE. 
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Y B. (GAMMA - 1.0) * UE/(S.O * GAMMA) 

DO 3 J c 1> MJMAX 
AMPL(J) B y 
PHASE(J) B 0.0 
3 CONTINUE 
GO TO 7 
B CONTINUE 

IF CNEDATA .EQ. 0) NZDATA « 5 
IF (NZDATA .EQ. 1) NZDATA « 7 
DO 6 I B 1« NJMAX 

READ (NZDATA, 5003) J# WIFLCJ), PHASE(U) 
e CONTINUE 

IF CNOZNLl .NE. 1) GO TO 7 

DO 710 1 B 1, NJMAX 

READ (NZDATA, 5003) J, GNOZ(J) 

710 CONTINUE 

7 DO 10 I B 1» NJMAX 

READ (5, 500a> J, L(J>, M(J), N(J), NS(J), NAME(J) 

10 CONTINUE 

DO IS J *= 1, NJKAX 

THETA *= PHASECJ) * PI/ 180*0 

YR B AMPLCJ) * COSC THETA) 

YI »= AMFLCJ) * SINC THETA) 

YNOZCJ) B CMPLX(YR,YI) 

IS CONTINUE 

ZE B 8.0 ♦ RLD 
CZE = CMPLX(ZE,0.0) 

C6AM = CMFLXC GAMMA, 0.0) 

CAX B CGAM 

IF (KDKOPS .EQ. 1) CAX •= CGAM + < 1.0,0. 0) 

l(ti)tljt+1^1<t^^c* + + 4 *+***^*****# + ** + **«t****#>|t*+***+:+* + ** + + **A************* 

ASSIGN ARRAYS FOR ROOTS OF BESSEL FUNCTIONS* 

DO 20 J B NJMAX 

IF C(M(J) .EG. 0) .AND. (N(J) .EQ. 0)) GO TO 15 
MM e M(J) + 1 
NN B N(J) 

S(J) B EJH00T(MK,NN) 

SJCJ) B RJVAL(MM,NN) 

GO TO S5 
15 S(J> B 0.0 
SJ(J> B l.O 
25 SSQ B S(J) * S(J) 

CSSQ(J) = CMFLXtSSQj 0.0) 

20 CONTINUE 

iJr + iJotoloC********#**^**^*******^*^**"^******************************** 

CALCULATE AXIAL ACOUSTIC EIGENVALUES. 

FIND MAXIMUM VALUES OF LCJ), K(J), AND N(J). 

KN B 0 
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LMAX = 0 
KMAX e= 0 
NI^AX c 0 

DO 30 0 = 1» MOM AX 


IF 

CLCd) 

• GT- 

LMAX 5 

LMAX 

S 

LCJ5 

IF 

CMCJ5 

.GT* 

MM AX 5 

MMAX 


MCJ5 

IF 

CNCd) 

.GT. 

NMAX5 

NKAX 

C; 

NCd5 

IF 

CNCJ) 

.NE. 

N< 1) ) 

KN = 

1 



30 CONTINUE 

LMAX = LKAX + I 
MMAX = MKAX + 1 

COMPUTE EIGENVALUES. 

DO AO J e 1, NJMAX 
LL c LCO) 

SMN •= S<d> 

YAMPL B AMPLCJ) 

YPHASE c PHASE<d> 

CALL EIGVALCLLj. SMN» GAMMA# ZE/ YAMPL* YPHASE* CEiLT> 

B(J> *= CESLT 
EC<J> « CONJGCCBSLT) 

AO CONTINUE 

CALCULATE LINEAF. COEFFICIENTS- 

CALCULATE THE NUMBER OF LINEAR COEFFICIENTS. 

KCOEFF = A 

IF fNOZKLj -EC. U NCOEFF = 5 
NCFMl *= NCOEFF- 1 

DO 100 NO = 1#’ NJMAX 
DO 100 NP = l> NdMAX 

ZERO COEFFICIENT ARRAYS. 

DO 105 KC *= I# NCOEFF 
cccKC#Nd*NP) *= CO. 0*0. o: 

105 CONTINUE 

ORTHOGONALITY PROPERTY OF TANGENTIAL EIGENFUNCTIONS- 

IF C NSCNP) -NE. NSCNd> ) GO TO 100 

IF (MCNPJ .NE.. NCNd>> GO TO JOO 

IF CMCNdJ .EQ. 0> GO TO 112 

AZ = PI 

60 TO ISO 

112 IF < NSCNdJ -EO. 15 GO TO 100 
^ = 2-0 * FI 

ORTHOGONALITY PROPERTY OF RADIAL EIGENFUNCTIONS. 

120 IF (NCNP5 *NE. NCNd5) 60 TO 100 
IF (S(NP)5 125* 122* 125 

125 SGM = M<Nd5 * MCNd> 

SSO = SCNP5 * SCNP5 
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SJSQ *= SJtNJ) * SJCNJ3 

RAD = CSSO - SeMJ * SJ£Q/<2«0 * SSO) 

GO TO JS7 
12S RAD =0-5 

CALCULATE AXIAL INTEGRALS. 

127 DO 130 WORT = A 

CALL AXIAL! (NOFT# NR-. NJ# UE» ZE# ZC0ME> CRSLTJ 
AX(NOPT) *= CRSLT 
130 CONTINUE 

EVALUATE FUNCTIONS AT NOZZLE END. 

ZEJ * CC0SH<C1*BC<NJ)*CZE) 

ZEFl •= CCOSH<CI*E<NF) + CZE) 

ZEP2 = Cl * B<NP5 # CSlNHf C1+B(NP)*CZE) 

CAZ = CMPLX<AZ>0.07 
CRAP *= CMFLX<RAp*0.0> 

COEFFICIENT OF THE SECOND DERIVATIVE OF ACP). 

CC<l#NJ>NF). B AX<n * CAZ ♦ CRAP 

COEFFICIENT OF AtP). 

CC<S,NJ^NF>^ = <CSSO(NF>*AXU> - AX<S3 + ZEP£*ZEJ> + CAZ * CRAP 

COEFFICIENT OF THE FIRST DERIVATIVE OF A< P> . 

CC<3»N0»NP3 «= <CAX*AX<3) + < g. 0» 0* 03*AX< AJ 
1 + CGAmYNOZ<NF)*ZEPl*ZEJ) ♦ CAZ ♦ CHAD 

COEFFICIENT OF THE RETARDED DERIVATIVE OF A(P>. 

CC<A^NJ>NP3 = CGAK * AX<3) * CAZ ♦ CHAD 

IF CNOZNLl .NE. 1) GO TO 100 

COEFFICIENT DUE TO NOZZLE NONLINEAfil TI ES. 

CESG = J- <GANHA-1> * UE/2. 

CC<5#NJjNP) *= UE * CESQ * GNOZ(KP) ♦ ZEJ ♦ CAZ ♦ CHAD 
100 CONTINUE 

NORMALIZE LINEAR COEFFICIENTS. 

DO lAO NO = l» NUHAX 
CNOKM<NJ) « CC(UNJ>NJ) 

DO lAO KP c 1. NJMAX 
DO lAO KC = Ij NCOEFF 

CC(KC.NJ>NP,) = CCCKC.NU.NP) /CNORMCNJ) 
lAO CONIINUE 

♦ 3fi*it3t:*** + ** *♦»*#**#*♦<=**♦*****+ + ********♦**♦**« ******* 

COMPUTE NONLINEAR COEFFICIENTS* 

IF (KONLIN .EG. 0) GO TO A02 

Cl »= tCGAM • <1.0/0.0>) * <0*5.^0*0> 
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COMPUTATIONS OF BESSEL INTEGRALS VHEM ALL SERIES TERMS HAVE THE 
SAME, RADIAL MODE NUMBER N(J)* 

IF <KN .EO* 1) GO TO 170 
DO 150 MP = MM AX 

DO 150 MS If MMAX 
DO 150 MO »= If MMAX 
BESICKFj>MQ^MJ) = 0.0 
BESBCMP^MQ^MJJ = 0.0 
BES3CMP>MO»MJ3 = 0«0 
LI «= MF - 1 
L8 = MQ - 1 
L3 = MJ - 1 
LM = LI + LS 
LN = LI + L3 

MN = LS + L3 . „ „ . 

IF C<L3.E0.LM> »0R* <L£»E0.LN5 .OR* CL1*E0.MN)> GO TO 160 

60 TO 150 

160 IF (NWAX *EQ. 0) GO TO 165 
Ai = ROROOT<MP^ NMAXl ' 

A2 = ROROOTfMO^NMAXJ 
A3 “ RJH00T<MJ>NMAX> 

GO TO 167 
165 Al c 0«0 
A2 = 0.0 
A3 = 0.0 

167 CALL RADIAL< 1»L1>LS»L3j>A1>A2^A3^ RESULT) 

BESICMF^MQ^KJ) = RESULT 

CALL RADIAL<8*LI>LEjL3> A1^ A8j A3» RESULT) 

BES2(«P*M0»Mc)> RESULT 

CALL RADIAL< 3» L ULO^LS# Al> A8> A3> RESULT) 

BESSCMF^MG^KJ) = RESULT 
150 CONTINUE 

170 DO 200 NJ = NdWAX 
DO 200 NF = 1> NJMAX 
DO 200 NO = If NJMAX 

CDl CNJ^NP^NQ) = CO. 0*0.0) 

CD8<KJ*NP*NG) = CO. 0*0*0) 

CD3CNJ*NP*NC) - CO. 0*0*0) 

CDACNJ#NF*NO) = <C»0*0.0) 

DO SIO J = 1* S 

CALL AZIMTLCJ*NP*NO*NJ* RESULT) 

AZICJ) = RESULT 
TANINTCJ) = CKFLXC RESULT* 0* 0) 

810 CONTINUE 

IF CASrCi)) 280* 225* 820 
225 IF (A2K2)) 220* 200* 280 

820 IF CKN *EQ. 0) 60 TO 282 

LI = MCNP) 

L2 = KCNG) 

L3 = MCNO) 
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A1 = S<nP> 

AS e SCNQ) 

A3 B SCNJ} 

GO TO SAA 

ass MP a MCNP) + 1 

MQ = KCNQ> + 1 

MJ = M(NJ> + 1 

BADINT(I> = CMFLXCBES 1 CMF>MQ^M 03 ^ 0 « 0 > 
BADINTCSJ *= CMPLXCBES 2 C«F#MQ#MJ)- 0-03 
BAD 1 NTC 3 > *= CMPLXCBES 3 C«F*MQ>MJ 3 * 0 . 0 ) 

SAfl DO SAO J = l> 3 

IF CKN .EQ. O) GO TO 2 A 2 

CALL RADIAL CJj>Ll#L 8 /L 3 » Al. A 2 »A 3 # RESULT) 
RADINTCJ) e CMFLXC RESULT. 0 . 0 ) 

242 DO 240 NC c 1,4 

CALL AXIAL 2 CO.NC.NP.NQ.NJ.EE. CRSLT) 
AXINTCNC.J) «= CP.SLT 
840 CONTINUE 


DO 250 J *= 1.4 

T1 = G1 ♦ CSSG<NP) * AXINT<0. n 
T9 *= G1 * AXINTCJ.3) 

D1 «= AXINTtJ. 1) ♦ TAKINTt 1 ) + RAUNT<33 

D2 «= AXINT<J. 1 ) + TAN1NT<2) + RAD1NT<25 

D3 = AXlNTtU.e) * TANINT<1) * RADlNTt 1) 

D4 = CT2 - Tl) 'K TANINTCl) .* RADINT<1) 

DCOEF «= <0*5.0«0> * <D1 + D2 + D3 + D4)/CN0RM<NJ) 
IF <J .EQ. 1) CDICNJ.NF.NG) » (1.0.- 1*0) ♦ DCOEF 

IF <d -EQ. 2) CE2<Nd.NF.NG) - (1.0. 1.0) * DCOEF 

IF <J .EG. 3) CD3(NJ.NP.NG) .= <1.0. 1.0> + DCOEF 

IF <J .EQ. 4) CD4<Nd.KP.N0) *= (l.O.-l.O) * DCOEF 

250 CONTINUE 
200 CONTINUE 

CALCULATE COEFFICIENTS FOR EQUIVALENT REAL STSTEH. 

402 DO 350 NJ = 1. NUN AX 
NEVd c <2 ♦ Kd) 1 
NEWdl = KEVd + 1 
DO 350 NP = 1. NdMAX 
REVP ~ (.2 * NP) -- 1 
NEVPl = KEVF +. 1 

COEFFICIENTS OF LINEAR TERMS. 

CCR « REALCCCC l.Nd.NP) ) 

CCI = AIKAGCCCC l.NJ.NP) ) 

CHKEWd.NEWP) n CCR 
Cl (NEVd.NEVPl) = -CCI 
ClCNEWJl.KEVP) «= CCI 
Cl(NEWJl.NE\vPl) « CCR 
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CCR »= RE^l.CCCCKC+ UN0,NF3 ) 

CCI = AIMA6<CC<KC+ 1#NJ^NF) ) 

cckcjKEwj>nev:f> = ccf 

CCKC>NEWJ»NE.VP1 ) - -CCI 
CCKCjNEV/JI^NEVjF) t= CCI 
CCKCjNLVJI^NLVFI ) = CCR 
360 CONTINUE 

COEFEICIENTS OF NONLINEAR TER^S* 

IF (NONLIN .EC. 03 6C TO 350 

DO 370 NO = NJMAX 

NEWQ = C2 * NO - I 

NEWei *= NEWO +1 ' ■ 

CD IB = BEAL<Ct'l<NJ>NF^KQ3 ) 

CEII = AINAG<CmCNJ>NF*N03 ) 

CDSfi «= fiEALCCE2(NJ>NF>KO 3 
CE2I c AIMAG<CD2<NJ^NF>NG) ) 

CD3R = EEAL<CE3CKJ>NP/NO) 

CE31 *= A1NAG<CE-3(NJ^NF#NG>> 

CD4E = EEAL(CEA(Nd>NP.NG) ) 

CD4I = AIMA6(CBA<NJ>NE*K0) 

E<KEWd#NEfePi KEVO> = CDIB + CDSE + CE3R + GEAR 
D<NEWj NEVF>NEV&1 3 = -CEII + CESI - CE31 + CE4I 

E<KEViJ# NE.V'PUKEV03 = -CDII - CESI + CD3I + CE41 

E<NECJ>NEUF1^KEVG13 = -CEIE + CE2E + CE3E - CE4R 
ECNEV,'Jl»NEfeP>N£WO> *= CBJI + CESI + CE3I + CDAI 
EtNEk'dl^ NEWP^NEWQ1> = CDIR - CESR + CE3R - CEAP. 

EtNEV/dUNEWPi^NEkO) = CEIR + CESR - CE3E - CEAR 

D<NEWdUNEfcPl>NEVGl3 -CEII + CESI + CD3I - CDAI 
370 CONTINUE 
350 CONTINUE 

+ #**** + ^c+**=*:it: + + =J=>t'+ + + + =i'*++*** + + *** + +4+=K^= + ** + + **=)‘^**:*=*'=*=***+++ + ^***** 

CONFUTE COEFFICIENTS FOR THE EQUATIONS VHICR ARE DECOLFLEE 
IN THE SECONE DERIVATIVES. 

EO 405 KC = 1> NCOEFF 
KKAXCKC3 « 0 
405 CONTINUE 

CALCULATE INVERSE OF THE MATRIX ClCl#d>. 

ON AX «= NdMAX 
NdMAX = S * NJMAX 

V(13 == 1 

CALL GdR(Cl^MX8*MX2^NdMAX#0* SSOO^dCj. V3 

USE INVERSE TO CALCULATE DECOUPLED COEFFICIENTS. 

DO 410 NP t: 1, NJMAX 

LINEAR COEFFICIENTS. 

DO 420 Nd *= I* NdMAX 
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PO 420 KC * NCFMl 
TS(KCiNJ) = 0«0 
DO A20 K » I# NiJFiAX 
tS<KC*NJ) = TStKC/NJ) 

/teO CONTINUE 

DO i330 NO = 1/ NJKAX 
DO A25 KC *: 1, 3 
CCKC»NJj-NF) = TSCKC^NJ) 

/iBSVAL *= AESCC(KC>NJ>NP>> 

IF CABSUAL •GE« SMI) KMAXCKO 
425 CONTINUE 

IF <N0ZNL1 -NEo 1) GO TO 430 
C<4»N«J#NP) *: TSCA^NJ) 

ABSUAL c ABS<C<4#NJ^NP>) 

IF <ABSVAL .GE. SM3) KKAXC4) 
430 CONTINUE 


+ CHNJ*K) * C<KC#K/NP> 


KMAXtKC) 


KMAXC4) ♦ I 


C 

c 


440 


445 

415 


NONLINEAR COEFFICIENTS* 

IF (NONLIN «E6* 0> GO TO 410 
DO 415 NO » 1# NJMAX 
DO 440 NO « !-» NOKAX 
TSQtNJ) * 0*0 
DO 440 K * 1* NOMAX 

TSQ(N0> *: TSQ<NO> + CUNOjK) * D<K>NP^NO) 
CONTINUE 

DO 445 NO e: 1> NOMAX 
DCNOiNpjNQ) c TSO<NO> 

ABSUAL = ABS<DCNJ»NP»NQ)> ^ 

IF CABSUAL »GE* SM2) KMAX(NCOEFF) *= KMAXCNCOEFF) 

CONTINUE 

CONTINUE 


410 CONTINUE 

♦ ♦♦♦**♦«*****»*♦♦>(«****»•«♦***+>•=****+♦*+*♦******♦******************** 
OUTPUT* 


IF CNOUT *GE* 2) GO TO 455 


PRINTED OUTPUT 

VRITE (6»600l) CTITLECD# I = 1# 72) 

VRITE <6^6002) GAMMA* UE»RLD*ZCOWB 
IF CNDROPS *E0* 0) VRITE (6*6020) 

IF (KDROFS *£0* D VRITE (6*6021) 

IF (NOZZLE .B.0* 0) VRITE (6*6012) 

IF (NOZNLl .EQ* 1) GO TO 760 
VRITE (6*6028) 

VRITE (6*6004) 

DO 310 0 c 1* OMAX 

VRITE (6*6003) NAME(O)* 0* LCO)* M(0)* N(0)* 
I ■ S(0)> £0(0)* ECO)* TNOZCO) 

310 CONTINUE 
GO TO 765 
760 CONTINUE 

VRITE (6*6023) 


NSCO)* 


89 



o o o 


WRITE (6>6025> 

SrITE°C 6 » 6026 /*^SwE<J>- J- LCJ>' NtJJ* NSCd># 

^WRITE C6,6026J YN0ECJ3. GiiOZiJ^ 

770 CONTINUE 
765 CONTINUE 

IF CNONLIN •EG* 0> WRITE (6>60135 

OUTPUT OF LINEAR COEFFICIENTS* 

BO 320 KC = 1>. KCEK 1 

If <KC *Ee* 15 WRITE <6 j6005> 

IF tKC .EG* 25 WRITE C6 j6006) 

IF CKC *EQ. 35 WRITE <6^6007) 

IF CKC *E0* A5 WRITE <6^602A> 

WRITE <6*6008 5 < d* d = 1* NJKAX5 

WRITE <6*60145 

EO 320 Nd = 1* NJKAX i, 

WRITE <6*6009 5 Nd* < C<KC*Nd*NP>* NP = I* NdKAK) 

3S0 CONTINUE 

; OUTPUT OF NONLINEAR COEFFICIENTS* 

IF <NOKLIN *EG* 05 GO TO 452 
DO 400 Nd = 1* NdMAX 
WRITE <6*60105 Nd 

WRITE <6*60115 Cd* d s= 1* NdMAX) 

WRITE <6*60155 
BO 400 NF = 1* NdMAX 

WRITE <6*6009 5 NF* < D(Nd*NP* N05* NG = 1* NdMAX> 

/jOO CONTINUE 

A52 if <N0UT *E0* 05 60 TO 4 

C 

A55 IF <NOUT »EQ* 35 GO TO 480 

WRl TE' COEFFICIENTS ON FASTRANB^ FILL. 

WRITE <9*70015 GAKEOA* UE* ZE* ZCOMB* NEROPS* NdMAX* NOZNL 1 

^ITE°<9*70025 d* LCd5* M<J>* NCd5* NS<d5* SCd5* SdCJ5* 
j NAKE<J5 

A50 CONTINUE 
C 

BO 457 J = 1* dMAX 

WRITE <9*70065 d* YNOZCd)* B<d5 

A 57 CONTINUE 

IF (NOZNL 1 -NE* 15 GO TO 720 
EO 730 J = 1* UMAX 
write <9*70075 d* GN0Z<d5 
7 30 COfJTIKUE 
720 CONTINUE 
C 

EO 460 KC = 1^ 3 
WRITE <9*7003) KMAX(KC) 

DO 460 Nd I* NdMAX 
DO 460 NP = 1* NdMAX 
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ABSVAL *= ABStCCKC»NJjNPJ> ^ 

IF (ABSVAL .GE. SM1> VRITE <9*70045 NJ*NP* CCKC>NJ#NF5 

a60 continue 

IF <NOZNL.l .NE* 15 GO TO 464 
WRITE <9*70035 KKAX<45 
DO 462 NJ “ I* NJWAX 
DO 462 NF s 1* NJMAX 
ABSVAL *= ABS<C<4*NJ*NP)5 

IF CAESVAL -GE. SM35 WRITE <9*70045 NO* NP* C<4*NU*NF> 

462 CONTINUE 
464 CONTINUE 

WRITE <9*70035 KMAX<NCOEfF) 

IF CNONLIN .EG. 05 60 TO 4 

DO 470 NJ e 1* NJl<'AX 
DO 470 NP *: 1* NUMAX 
DO 470 NO «= 1, NJNAX 
ABSVAL = ABSCD<NJ*NP*NQ5) 

IF CAESVAL .GE. SM25 WRITE <9*?0055NJ» NF* NO* D(NJ>NF*Ne) 
470 CONTINUE 
GO TO 4 


PVNCRED CARD OUTPUT. 


460 PUNCH 7001 GAKNA* UE* EE* 

DO 482 0 c 1* UMAX 
PUNCH 7002 J* LCJ5* MCJ5* 
I NAHE<J> 

462 CONTINUE 


ZCOHB* NDROFS* NJMAX* NCZNLl 
N<J5* NSCJ5* S<J5* SJCJJ* 


DO 464 0 *= I* UMAX 
PUNCH 7006 J* YN0ZCJ5* E<J5 
484 CONTINUE 

IF (NOZNLl .NE. 1) GO TO 740 
CO 750 J *= I* UMAX 
PUNCH 7007 J*GN0Z<«J5 
750 CONTINUE 
740 CONTINUE 

DO 486 KC = 1* 3 
■ PUNCH 7003 KKAXCKC5 
DO 486 Nd = 1* NdKAX 
DO 466 NP *= 1* NdMAX 
ABSVAL » ABS<C(KC*Nd*NP5) 

IF C ABSVAL .GE. SK 1 5 PUNCH 7004 Nd* NP* 
486 CONTINUE 

IF <N0ZNL1 ®NE. 1) GO TO 490 
PUNCH 7003 KHAX<4> 

DO 492 Nd = I* KJKAX 

DO 49 2 NP = 1* NJMAX 

ABSVAL « ABS<C<4*Nd*NF>5 

IF < ABSVAL .GE. SM3> PUNCH 7004 Nd* NP* 


C<KC*NJ*NP) 


CC4*NJ*NP> 
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492 CONTINUE 
49 0 CONTINUE 
C 

PUNCH 7003 KKAX CNC0EFF3 
IF <N0NLIN .E6. 05 GO TO 4 
DO 488 NO *= 1j NJKAX 

DO 488 NP = NdMAX 

DO 488 NO = NJKAX 

ABSVAD = ABS<DCNJ^NP^NG) ) 

IF (ABSVAL *6&. EK25 PUNCH 7005 NJ> NF, NO.* D<NJ>NP^N05 
488 CONTINUE 
CO TO 4 

ERROB EXIT 

500 IF S!0» 510> 520 

510 JCC15 = ABS<dC( 1)5 
WRITE <6 j 6017) JCCl 5 
GO TO 4 

5S0 WRITE C6>6018) JCCl 5 
GO TO 4 
600 CONTINUE 

WRITE (6^60975 




FORNAT 

5000 FORNAT 

5001 FORMAT 

5002 FORMAT 

5003 FORMAT 

5004 FORMAT 

5005 FORM'AT 

6001 FORMAT 

6002 FORMAT 
1 

6003 FORMAT 

6004 FORMAT 
1 

6005 FORMAT 

6006 FORMAT 

1 

6007 FORMAT 
1 

6008 FORMAT 

6009 FORMAT 

6010 FORMAT 
1 

6011 FORMAT 

6012 FORMAT 

6013 FORMAT 

6014 FORMAT 

6015 FORMAT 

6017 FORMAT 

6018 FOF'MAT 
6020 FORMAT 


SPECIFICATIONS* 

(72A15 

(4F10.0j.2I 55 
( 51 5» IXj A4) 

(15# 2F 10.05 

(6155 

<2F10.0> 

< IHl# IX# 72Ai//'5 

C2X#8HGAM.MA = # F 5. 8# 5X# 5HUE = #F5.2# 5X> 6HL/D ~ #F8i5# 
5X>8HZC0MH *= »F5.2/5 
(2X#A4>5I 5# 4F10.5#2F11. 5/5 

(2X////EX# 29HNAKE J L M N NS# 7X# 3H£MN> 3X# 
7HJM( SMN5# 7X# 3HEFS# 7X# 3HETA#8X# 2HYR#.8X# 8HYI //5 
<IHI#45H DECOUPLED COEFFICIENT OF B(P5* C< 1# J# P)///) 

CIH1#44R DECOUPLED COEFFICIEin’ OF THE DERIWATrUE OF# 

6H B(P):# 5X#8HC(2#.J#P)///5 
(1HI#39H DECOUPLED COEFFICIENT OF THE RETARDED# 

BOH DERIVATIVE OF B( P) T # 5X# 8HC( 3# J# P) /// ) 

<7X# IHP# I8#9I 125 
(8X//8X#I3#3X# 10F1&.65 

(lHi#4SH DECOUPLED COEFFICIENT OF ECF5 * DBC05/DT# 

19H IN EQUATION FOR B( # I S# IH 5 /// 5 
(7X# IHG# I8>91 185 
( 8X# 19HQUAS1- steady NOZZLE/) 

<2X//8X#84hLlNEAR COEFFICIENTS ONLY) 

(4X# 1HJ5 
(4X#1HF) 

(1H1#31H OVEFFLOV DETECTED# LAST ROV = # I 5) 

(1H1#34H singularity DETECTED# LAST ROW *= # 1 55 
(2X# ’'DROPLET MCMENTW SOURCE NE&LECTEL’V) 
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60S1 

6022 

6023 

6024 

6025 


6026 

6027 

7001 

7002 

7003 

7004 

7005 

7006 

7007 


FOEt'AT C2X* "DROPLET KOKENTDM SOURCE INCLUDED"/) 

FORMAT C2X, "NOZZLE NON-LINEARITIES NEGLECTED"/) 

FORMAT t2Xi "NOZZLE KGNLINEARITI ES INCLUDED"/) 

FORMAT ORl^" DECOUPLED COEFFICIENT DUE TO NOZZLE"/ 

" NONLINEARITIESj"/ 5X/8hCC4/J/P)////) 

FORMAT C2X////2X, 29HNAKE • J L « N NS# 7X> 3HSMN/ 3X> 
7hdM < SMN ) # 7X # 3H EPS# 7X# 3HETA# 8X# 2KV R# 8X# 2H7 I # 

S 8X#2HGR#8X# 2R6I-//> 

FORMAT C2X#A4#5I5#4F10.5#4F11.5/) 

FORMAT ClHl) 

FORMAT C4F10«5#31S) 

FORMAT C5I5#2F10.£# iX#A4) 

FORMAT CIS) 

FORMAT (2I5#FJ5*6) 

FOFsMAT (3I5#F15*65 
FORMAT <15#4F10.5) 

FORMAT <I5#2F10.5) 

END 
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SUBROUTINE EI6VAL<L^ ShN# GAMMA^Z E*YAKFL>T FHASEj RESULT) 

COMPLEX RESULT • ' ■ ’ , 

COMMON /ELKl/ GSB# AESO# ALPET- SMNS0 

THIS SUEROUTINE COMPUTES THE COMPLEX AXIAL ACOUSTIC EIGENVALUES 
FOR A CYLINDRICAL CHAMBER WITH A NOZZLE AND STORES THEM IN 
RESULT. 

THE EIGENVALUES ARE COMPUTED BY MEANS OF NEWTONS METHOD* 

THE INPUT PARAMETERS APE AS FOLLOWS 
LIS THE AXIAL MODE NUMBER. 

SMN IS THE DIMENSIONLESS ACOUSTIC FREGUENCY. 

GAMMA IS THE SPECIFIC HEAT RATIO. 

ZE IS THE LENGTH- TO-RADI US RATIO. 

YAMFL IS THE NOZZLE AMPLITUDE FACTOR. 

YFHASE IS THE NOZZLE PHASE SHIFT IN DEGREES- 

PI « 3.1A15927 
ERR = 0-0000001 

IF (YAMFL) 5*60»5 
CALCULATE CONSTANTS. 

5 PHASE = YFHASE * FI/180.0 
ALPHA = YAMFL + CO S( PHASE) 

BETA = YAM PL ▼ SIN (PHASE) 

6S0 = GAMMA * GA>:HA 

ABSe = (ALPHA * ALPHA) - (BETA * BETA) 

ALBET ALPHA * BETA 
SMNSQ *: SMN * SMN 

ASSIGN INITIAL GUESS FOE EIGENVALUE. 

IF <L ,.EG. 0) GO TO AS 
RL «= L 

PHI = Pl/e-0 + PHASE 
XK *= RL PI/ZE 
A «= YAMFL/ZE 
XO = XM + A+COSCPHI) 

YO s AX'SINtPHl ) 

60 TO A7 
AS CONTINUE 

YPHl = YPHASE 

IF (YFHASE -GT. 180) YFHI = YFHASE - 180. 

IF (YPHASE .LT. 0) YFHI = YPHASE + 180- 

IF (YAMPL .LT. 0-1) 60 TO UO 

IF (YAMFL .LT. O.A) GO TO 120 

IF- (YAMPL .LT. 0.6) 60 TO 150 

IF (YAMPL .LT. l.S> GO TO 160 

XO = 1.0 + YAMPL 

GO TO 170 

160 XO « 1 .as * YAMPL 
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170 IK CYFHI -LEo 30*5 TANFSI = -0«4 

IF <YPHl*67«30« .AND* YFHI.LE.60.3 TANPSI = -0«2 
IF. <YFHI.GT*60« .AND* ,YFHI .LE. 180»5 TANPSI «* 0*0 
IF CYPHI *GT» l£0« *AKD* YFBl .LE* 1 50* I TANPSI » 0*8 
IF <YPHI .GT. 150.) TANPSI ** 0«A 
GO TO 140 

150 XO “ 2*0 * YAMFL 

IF (YFPI .LE. 30.) TANPSI » -0*6 

IF (YPH1*6T*30* .AND* YFHI *EE* 60* ) TANPSI *= -0*3 
IF CYFHI »CT*60* .AND* YPHl »LE« 120* > TANPSI ~ 0*0 
IF CYFHI .GT. 120* .AND. YPHl *LE. 1 50* ) TANPSI ^ 0*3 
IF CYFHI *GT* ISO*) TANPSI » 0*6 
GO TO 140 

no XO *= 5* \ YANFL 
60 TO 130* 

120 XO = 3* * YAMPL 
130 CONTINUE 

IF CYPHI *LE. 30.) TANPSI = -0*75 
IF CYPHl .GT.30* *AND» YPHI.LE.60*) TANPSI « **0»4 
IF CYFHI*6T«60* «AND* YPHl *LE* 1 80* > TANPSI *= 0*0 
IF CYPHl *GT* 180* *AND* Y PHI *LE* 1 SO* ) TANPSI ** 0*4 
IF CYPHl *GT*. 150.) TANPSI « 0*75 
140 CONTINUE 

YO - XO + TANPSI 


ITEBATION USING NE^^TONS METHOD FOE A SYSTEM OF TfcO EQUATIONS 
IN TV'O UNKNOWNS* 

47 LI = 0 
X = XU 

Y = YO 

40 CALL FCNSCX>YjZEjFjG.» F^>FY.fGX..6Y) 

IF CLl *E0* 40) GO TO 50 
RJFG » CFX + GY) - CGX ♦ FY) 

IF CRJFG) 20* 30* 20 
20 DEL TAX cC-F*GY + G* FY)/EJFG 
DEI-TAY = C-G * FX + F * 6X)/B*JF6 
LI = LI + 1 
X = X + DELTAX 

Y = Y + DELTAY 


C 

c 


c 

c 


c 

c 


TEST FOR CONVERGENCE. 

IF CABSCDELTAJO *GE. ERE 
GO TO 10 


.OB. ABSC DELTAY) 


.GE* EBR) 


WARNING MESSAGES 
30 WRITE C6>6005) 

GO TO 10 

50 WRITE <6*60065 
-GC TO 10 

CASE OF HARD WALL CYAMPL = 05* 
60 EL = L 


GO TO 40 
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X « KL ♦ PI/ZE 
Y = 0»0 

10 RESULT = CMFLXCX»Y) 

FORMAT SPECIFICATIONS. 

6005 FORMAT < £X//2X* 16HUAC0EI AN IS ZERO//> 

6006 FORMAT < 2X//2X^ SSKFAILEP TO CONVERGE IN AO ITERATIONS//) 
RETURN 

END 
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SUBHOUTI N E FCN S C Z FX^ FY^ GX j- GY > 

C 

C THIS SUBROUTINE COMPUTES THE FUNCTIONS FCX,Y> AND G<X>Y) 

C AND THEIH PARTIAL DERIVATIVES WITH RESPECT TO X AND Y* 

C 

COMMON /BLKl/ GSQ> ABSQj ALBETj SMNSQ 
C 

C COMPUTE THE TRIGONOMETRIC FUNCTIONS^ THE HYPERBOLIC FUNCTIONS 

C AND THEIR SQUARES* 

C 

1=1 

ABGX = ZE 4= X 
ARGY = ZE * Y 
10 SX = SINCARGX5 
CX = COSCARGX) 

SHY = SINK<ARGY3 
CHY = CQSHtARGYJ 
IF <I *EQ. 2) GO TO 20 
SXSQ = SX 4 SX 
CXSQ = CX * CX 
SHYSQ = SHY + SHY 
CHYSQ = CHY * CHY 
AEGX = 2.0 4= ARGX 
ARGY = 2*0 4= ARGY 
,1 = 2 
GO TO 10 
C 

C COMPUTE TRANSCENDENTAL FUNCTIONS AND THEIR DERIVATIVES 

C 

20 FF = tsxse * CHYSQ) - <CXSQ ^ SHYSQ) 

GG = CCXSQ * CHYSQ) - C SXSQ * SHYSQ) 

HH = 0*25 4= SX * SHY 
FFX = ZE * SX * CHY 
GGY = ZE * CX * SHY 
FFY = -GGY 
GGX = -FFX 
HHX = 0*5 * GGY 
HHY = 0*5 * FFX 
C 

C COMPUTE FACTORS 

XYSQ = <X * X)' - <Y 4= Y) 

XY = X 4= Y - 
SMNXY = SMNSQ + XYSQ 

FI = (ABSQ * SMNXY) - (4*0 * ALBET * XY) 

F2 = (ALBET 4= SMNXY) + (ABSQ * XY) 

G: = (ABSQ * SMNXY) + (4*0 4= ALBET * XY) 

FXl = (2.0 * X 4 = ABSQ) - C4.0 4: ALBET * Y) 

FX2 = (2.0 * X 4= ALBET) + (ABSQ * Y) 

FYl = (-2.0 4 = y 4: ABSQ) - <4.0 * ALBET +- X) 

FY2 = (-2.0 4= Y 4= ALBET) + (ABSQ * X) 
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GXl - C2.0 * X * ABSQ) + (4-0 * ALBET *' Y) 
GYl = C-2.0 * y * ABSQ) + C4.0 * ALBET * X> 

COMPUTE F<XjY) and GCXiY) 

F'= (XYSQ * FF> - (4.0 * XY =* HH) 

1 + GSQ * <<F1 * GGJ + (4.0 * F2 * KH)) 

G = (XYSQ * HK) + (XY * FF> 

1 + GSQ * C(F2 * GG) - CGI * HH)3 

COMPUTE THE PARTIAL DERIVATIVES OF F AND G 

FX = C2.0 =f= X * FF) + CXYSQ * FFX) 

1 -4.0 + ((Y * HH) + CXY * HHX)) 

2 + GSQ * CCFXl + GGJ + CFi * Q6X3 

3 + C-4.0 * FX3 * HH) + <4.0 * F2 ♦ HHX) > 

FY = C-2.0 4: y 4! FF) + CXYSQ + FFY) 

1 -4.0 * (CX * HK) + CXY *■ HHY)) 

2 + GSQ * CCFYl * GG) + CFI * GGY) 

3 + C4.0 * FY2 HH) + <4*0 =t= F2 * HHY)) 

GX = C2.0 * X HH) + (XYSQ =f= HHX) 

1 + CY * FF) + (XY * FFX3 \ 

2 + GSQ * (CFX9 =i' GG) + CF2 * GGX) 

3 -C6X1 * HH) - (Gl * HHX)) 

GY = C-2.0 + Y * HH) + CXYSQ =i= HHY) 

1 + CX =f= FF) + -CXY * FFY) 

2 + GSQ * CCFY2 + GG) + CF2 * GGY) 

3 -CGYl * HH) - CGI * HHY))' 

RETURN 

END 
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SUBROUTINE AXI/iH <NOPT>NP#NJ^ UE*2E;EC0MB* RESULT) 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


n 


THIS SUBROUTINE CALCULATES THE INTEGRAL OVER THE INTERVAL 
<0>2E) OF THE FOLLOWING FUNCTIONS ACCOHEIN6 TO THE VALUE 
OF NOPT 

NOPT «* 1 Z(NP) * ZC(Nd) 

NOFT * R ZPP<NF> * ZCCNd) 

NOPT =3 UP * ZCNP) * ZCCNd) 

NOPT «= fl U + ZP<KF> * ZCCNd) 

IN THE ABOVE ESUATIONSS 

ZCNF) IS THE axial ACOUSTIC EIGEWFUNCTI OK OF INDEX NF* 

Z<KJ> IS TKF AXIAL ACOUSTIC EIGENFUNCTION OF INDEX KU< 

ZC IS THE COMPLEX CONJUGATE OF THE AXIAL EIGENFUNCTION. 

ZP AND ZFF ARE THE FIRST AND SECOND DERIVATIVES OF THE 
AXIAL EIGENFUNCTIONS RESPECTIVELY, 

0 IS THE STEADY STATE VELOCITY DISTRIBUTION AND UP IS ITS 
AXIAL DERIVATIVE. 

THE VELOCITY DISTRIBUTION IS COMPUTED BY THE SUBROUTINE UBAR. 


parameter MX = 5 
REAL MAG 

COMPLEX CI» CZE, BP* BJ* Tl* T2f CH* FI* F2* F3* CZ* ARG* 
I SI* £2* S3* RESULT* FUNCTC 500) * BCMX) 

COMMON B 

Cl - CO* 0*1.0) 

C2E e CMPLXC2E*0.0> 

BP = B(NP) 

EJ !s CONJGIBCNJ)) 

IF CNOPT .GT. 2) GO TO 50 

CALCULATE INTEGRALS BY MEANS OF ANALYTICAL EXPRESSIONS FOR 
NOPT tr 1 AND NOPT » S. 

ARG = <BP + BJ) * Cl 
MAG = CABSCABG) 

IF <MA6) 20* 25* 20 
SO Tl = CSIKHCARG*CZE)/ARQ 
60 TO 30 
85 Tt <= CZE 

30 ARG «= CBP - BJ> ♦ Cl 
MAG = CABSCARG) 

IF <MAC) 35* <10* 35 
35 TS = CS1NHCARG*CZE)/ARG 
GO TO AS 
AO TS e CZE 

AS RESULT « <TI + T2) * CO* 5*0.0) 

IF CNOPT .EO. 2) RESULT s -BCNP> * ECNF) * RESULT 
GO TO 100 
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NUMERICAL EVALUATION OF INTEGRALS FOR NOPT = '3 AND' N OPT = 4* 

COMPUTE STEP SIZE FOR SIMPSON INTEGRATION* 

50 N = 50 ■ ' 

RN = N 

RESULT = <0*0>0*0) 

IC « ZCOMB 
IC = 2 - IC 

DO 90 J c 1, IC 

IF CJ «EQ* n H = ZCOMB * ZE/RN 

IF *EQ. 2) H = (1.0 - ZC0MB3 * ZE/RN 

IF CJ .EQ. 1) ZO = 0.0 

IF CO .EQ. 2) ZO = ZCOMB ^ ZE 

NPl = N^,+ 1 

CH = CMPLXCH^O.O) 

COMPUTE INTEGRANDS. 

DO 60 I « 1> NPl 
STEP = I “ _1 - 
Z = (STEP * H3 + ZO 

IF CCI.EQ.n .AND. (J.EQ.SJ5 Z >= Z + H/100.0 
IF CNOPT *EQ. 33 CALL UBARC 2> UE^ZE»ZCOMB^ Z>F) 

IF CNOPT -EQ- A3 CALL UBARC 1> UE^ZEjZ COMB#Z> F3 
FI = CMPLXCF>0.03 

CZ == CMPLXCZ^0.03 • : 

ARG = Cl * BP 

IF CNOPT .EQ- 33 F2 = CCOSHC ARG*CZ 3 
IF CNOPT .EQ. 43 F2 = ARG * CSINH(ARG=«=CZ3 
ARG = Cl + Bd 
F3 = CC0SHCARG*CZ3 
FUNCTCI3 = FI * F2 * F3 
60 CONTINUE 

PERFOB^J SIMPSON INTEGRATION. 

NMl s N - I 

51 = FUNCTCI3 + FUNCTCNP13 

52 = CO. 0^0. 03 
£3 « CO. 0^0.0) 

DO 70 I = 2# N, 2 . 

52 = 52 + FUNCTCI3 
70 CONTINUE 

DO 80 I = 3^ NMl, 2 

53 - S3 + FUNCTCI3 
80 CONTINUE 

RESULT = RESULT + 

1 CH =*= CSl + (4.0*0.03=CS2 + C 2. 0> 0-03* S33 /C 3. 0^ 0. 03 

90 CONTINUE 

100 CONTINUE • 

RETURN 

END 
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SUBEOUTI NE AXIAL2CNOFT*N CON N N Q* N Z E» EESUL T ) 


C 

c 

C THIS SUBEOUTINE CALCULATES THE INTE6BAL OVER THE INTEEV/tt, 

C CO*ZE3 OF THE FOLLOWING FUNCTIONS ACCORDING TO THE VALUES 

C OF NOFT AND NCONJ 

C 

C FOB NCONJ «= 1 AND 

C NOPT « i ZCKF3 ♦ Z<NQ) ♦ ZC(NJ) 

C NOPT = 2 ZP<NP3 * ZF(NO) * ZCCNJ> 

C NOPT *= 3 ZFPtNP) * Z<NQ) * ZC(NJ> 

C 

C FOB NCONJ *= e AND 

C NOPT » 1 Z<NF3 + ZC<N03 * ZC(NJ) 

C NOPT *= fi ZFCNPJ * ZFCCNOJ * ZCCNJ3 

C NOPT «= 3 ZPP<NP> * ZCtNG) * ZCCNJ3 

C 

C FOE NCONJ *= 3 AND 

C NOPT * 1 2C<NP) * ZCNQ) * ZC<NJ> 

C NOPT = e ZPC<NF) * ZFCNG) * ZC(NJ) 

C NOPT = 3 ZFPCCNP) * ZCNO) * ZCtNJ) 

C 

C FOR NCONJ = A AND 

C NOPT « I ZC(NP) + ZC(NG) + ZCCNJ3 

C NOPT e S ZPC<NP> ♦ ZFC<NG3 * ZC<NJJ 

C NOPT e 3 ZFPCCNP) * ZCCNQ) ♦ ZCCNJ) 

C 

C IN THE ABOVE EQUATION St 

C Z(NP)i Z(NQ)> AND ZCNJ) ABE THE AXIAL ACOUSTIC EIGENFUNCTIONS 

C AND iiPf NO> AND NJ ARE IHEIR INDICES* 

C ZP IS THE FIRST DERIVATIVE OF THE AXIAL EIGENFUNCTIONS. 

C ZPP IS THE SECOND DERIVATIVE OF THE AXIAL EIGENFUNCTIONS* 

C ZC AND ZPC ARE COMPLEX CONJUGATES OF Z AND ZP HESFECTl VELY- 

C 

C 

PARAK.ETEB NX = 5 
REAL KAG 

COMPLEX Cl* PF* CZF* BF* BG* BJ* SUM* RESULT* 

1 ARGcii)* FtWCTCA)* BCMX) 

COMMON B 

CALCULATE INTEGRALS BY MEANS OF ANALYTICAL EXPRESSIONS. 

Cl « C0.0*l*0> 

CF B C0.25*0*0> 

CZE <= CMPLXCZE*0-0) 

BP = ECNP) 

BQ = BCNQ) 

BJ = CONJGCECNJ)) 

IF CCNCONJ.EQ.fi) .OR. CNCONJ* EG. A) ) BG *= CONJGCEG) 

IF CNCONJ .GT« S) BP « CONJGCBF) 

ARGCl) = CBP + BQ + BJ) * Cl 
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ARGC23 = (BP + BQ - BJ) ♦ Cl 

ARGC3> =. (BP - ;BQ + BJ) * Cl 

ARGC4) = (BP - ’bQ - BJ) * Cl 

DO 10 J = 1^4 
MAG = CABSCARG(J)> 

IF CKAG5 12^- 15j 12 
12 FUNCTCJD = CSINH<ARG<J1*CZE>/ARGC 
GO TO 10 

15 FUNCTCJ) = CZE 
10 COKTINUE 

IF CNOPT *EQ. 2D GO TO 30 

SUM = FUNCTCID + FUNCTC2D + FUNCTC3D + FUNCT<4) 
RESULT = CF + SUM 

IF CNOPT .EQ. 3D RESULT = -BP BP i RESULT 
GO TO 50 

30 SUM = FUNCTCID + FUNCTC2) - FUNCTC3D - FUNCTC4) 
RESULT •£= -CF * BP * BQ * SUM 
50 CONTINUE 
RETURN 
END 
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SUBROUTINE fiZIKTLCNOPT>NP»NQj NJj HESULT> 

C 

PARAMETER MX = 5 
DIMENSION NFCNC3>> S6C2> 

COMMON /BLK8/ MCMX)^ NS(MX> 

C 

C ^t^i^^Hfit:****************************************^******************* 

c 

C THIS SUBROUTINE CALCULATES THE INTEGRAL OVER THE INTERVAL 

C CO* S*PI) OF THE FOLLOWING FUNCTIONS ACCORDING TO THE VALUE 

C OF NOPT 

C 

C NOPT *= 1 THCNP> * THCNG) * THCNU) 

C 

C NOPT 2 ThPCNP) ♦ THPCNQ) ♦ TH<NJ> 

C 

C IK THE ABOVE EQUATIONS! 

C TH<NP>* THCNQJ* AND THCNd> ARE THE TANGENTIAL EIGENFUNCTIONS 

C AND NP* KQ* AND NO ARE THEIR INDICES. 

C THP IS THE DERIVATIVE OF THE TANGENTIAL EIGENFUNCTIONS- 

C 

C IF NS = 1 TH .= SINCM+THETA) 

C IP NS = 2 TH«= COSCK*THETA5 

C 

C 

RESULT = 0.0 
FACTOR *= 1.0 
PI a 3.1A159 27 


C 


C 


10 


DISTINGUISH BETWEEN 
DO 10 Kl' = 1* 3 
NFCNCHl) = 1 
CONTINUE 


SINES AND COSINES- 


IF CNSCNJ) .E&.2) 
IF (NOPT .EQ. 8> 
IF CNSCNP) .E£.2> 
IF CNSCNGJ .EQ.2> 
GO TO 30 

20 If CNStNP) .EQ. I ) 
IF (NSCNO.EQ. 1> 
DO AO Kl = 1*2 
SGCK1> = 1.0 
IF (NFCNCKp 

40 CONTINUE 

factor SG< I > * 


NFCNC3) = 2 
CO to 20 
NFCNC I) e a 
NFCK(2) = 2 

NFCNC n = 2 
NFCNC2) *= 2 


FG. 1) SGCKl) = -1.0 


SGC2> ♦ MCNP> + MCNG) 


30 NSUM - 0 

DO 50 Kl .= 1* 3 
NSIM «= NSUM + NFCN(Kl) 
50 CONTINUE 
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IF CCNSUM .EQ. 35 .OR. CNSUM .EQ. 555 60 TO 60 

IF <NSIW .EQ. 45 GO TO 70 ‘ ' 

IF CNSUM .EQ. 65 GO TO 60 

70 KOPT “ 2 

IF CNFCNC15' .EQ. 25 ' GO TO 72 
GO TO 74 
72 LL = MCNP5 
M« = MtNQ5' 

NN = MCNvJ5 
GO TO 90 

74 IF (NFCNC25 .EG. 25 GO TO 76 
GO TO 78 
76 LL = MCNQ5 
MM = MCNP5 
NN = MCNJ5 
GO TO 90 
78 LL = MCNJ5 
MM = MCNP5 
NN = MCNQ5 
GO TO 90 

60 KOPT - 1 
LL = MCNP5 
MM = MCN05 
NN i= MCNJ5 

COMPUTE VALUES OF THE INTEGRALS. 

90 IF {<LL.NE.05 -AND. CMM.KE.'o -AND. CNN-NE.055 ‘60 TO 101 
GO TO 103 

101 LM = LL + MM 

LN LL + NN ... 

MN = MM + NN / ^ 

IF {<NN.EQ.LM5 -OR. CMM.EQ.LN55 . RESULT = PI/2.0 
IF CLL .EG. MN) GO TO 102 
60 TO 104 

102 IF {KOPT .EQ. 15 RESULT = PI/2.0 
IF <KOPT -EO. 25 RESULT = '-FX/2.0 

GO TO 104 • ’ 

103 IF <{LL.EQ.05 .Al-JD. {MM.E0.05 .AND. CNN.EQ.05 > ’ GO TO, lOS 

IF ((KOPT. EQ. 15 .AND. (NN.EQ.05 .AND* CLL.EO.MM55' RESULT = PI 
IF ((K0PT.EQ.15 .AND. (MM.EQ.05 .AND. (LL.EQ.NN55 RESULT = PI 
IF <(LL .EQ. 05 .AND* (KM .EQ. NN55 RESULT = El. . 

GO TO 104 

105 IF (KOPT -EQ. 15 RESULT = 2.0 * FI 

104 CONTINUE 

RESULT = FACTOR * RESULT 
60 CONTINUE 
RETURN 
END 
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SUBROUTINE RADI ALCNOPTjL«MjNjO*Bj Cj RESULT? 

THIS SUBROUTINE CALCULATES THE INTEGRAL OVER THE INTERVAL 
C0>1> OF THE FOLLOWING PRODUCTS OF THREE BESSEL FUNCTIONS: 

NOPT == 1 JLCA+R? * JM<E+R? ^ JN<C*R? * R 

NOPT = 2 JLCA*R? * «JM<B*R) * JN<C*B?/R 

NOPT = 3 JPL<A*R? * UPMCB^R) * JNCC+R5 * R 

UL IS THE BESSEL FUNCTION OF FIRST KIND OF ORDER L 

JPL IS THE DERIVATIVE OF JL WITH RESPECT TO R 
L* N ARE NON-NEGATIVE INTEGERS 

A» Bt C ABE REAL NUMBERS 

DIMENSION FUNCT(200) 

DOUBLE PRECISION DN> DH# DSTEP> DH> ARG1> ARG2> AR63* 

1 BESIj BES2^ BES3^ BESH> BESL# PROD^ 

2 FUNCT# BESLIM# SI* SZ* S3 

NN 100 
DN s= NN 
DH = 1.0/DN 
NPl ■= NN + 1 

DO 10 I = NPl 
DSTEP =1-1 
DR = DH * DSTEP 
ABGl = A * DB 
ARG2 = B * DR 
ARG3 = C ={= DR 

CALL JBESCN>ARG3^BES3jS500? 

IF CNOPT «EQ» 3) GO TO lOJ 
CALL JBESCL,ARG1 j.BES1*$500? 

CALL JBESCM>ARG2^ BES2y$500? 

GO TO 102 

101 IF CL .EQ. 05 GO TO 103 

CALL JBESCL+1^ARG1^BESH>S500? 

CALL JBESCL-I^ ARGl^BESL^SSOO) 

BESl = A * <BESL - BESH)/2*0 
60 TO 104 

103 CALL JEESC 1^ ARG IjBESIj S 5005 
BESl = -BESl ^ A 

104 IF CM .EO. 05 GO TO 105 
CALL UBESCM+ 1^ ARG2j EESHj S5005 
CALL JBESCM-1^ARG2^EESLj S5005 
BES2 = B * CEESL - BESH5/2.0 
r;n ta mft 


105 



105 CALL JBES( ARGSjBES 2>S5003 
BESS = -BESS B 
102 PEOD = BESl * BESS * BESS 

IF ><NOPT *E0.’2> GO TO 110 
FUNCTdJ = PROD * DR 
GO TO 10 

110 IF <I .EQ. n GO TO 111 
FUNCTCI) = FROD/DR 

GO TO 10 

111 BESLIM =0.0 

IF CCL.EQ.ll .AND* CM*KQ.Q3 .i-uNu.- H t 

IF C<L.EQ.OT .AND.- CM.EQ.ll. .AND. CN.EQ.O)) BESLIM 

IF CCL.E0.03 .AND'. CW.EQ.OJ .AND. CN.EQ.U) BESLIM 

FUNCTUl = BESLIM 
10 CONTINUE 

NMl = NN - 1 

51 = FUNCTCI I. +.FUNCX(NP1 ) 

SS = 0*0 ' 

£3 = 0.0 

DO 20 I Js 2 ^ NN> 2 

52 = SS + FUNCTCI > 

SO CONTINUE 

DO 30 I = 3^ NMl, 2 

53 = S3 + FUNCTCI) 

30 CONTINUE 

RESULT = DH * (SI + 4.0*S2 + 2.0+S3)/3.0 
GO TO 501 

500 TOITE C6# 6000) 

6000 FORMAT C IHl# lOHERROR USES) 

501 CONTINUE 
RETURN 
END 


lo6 


A/S.O 

B/2.0 

!c>2.0 
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SUBROUTINE UBARCNOPT^ UK. ZE^ZCOMB^Z^ RESULT) 

THIS SUBROUTINE CALCULATES THE STEADY STATE VELOCITY 
DISTRIBUTION FOR L’NIFOEMLY DISTRIBUTED COMBUSTION COMH-ETED AT 
Z = ZCOKB * ZE WHERE: 

UE IS THE EXIT MACH NlSyiBER. 

ZE IS THE DIMENSIONLESS LENGTH. 

Z IS THE AXIAL COORDINATE. 

IF NOPT = 1 THE DISTRIBUTION IS CALCULATED. 

IF NOPT = 2 THE DERIVATIVE IS CALCULATED. 

IF NOPT = 3 THE SECOND DERIVATIVE IS CALCULATED. 


ECZ = ZCOMB * ZE 



GO 

TO 

< 10>20j30) 

j NOPT 


10 

IF 

<z 

• LE. 

ECZ) 

RESULT = 

UE * Z/EGZ 


IF 

cz 

-GT. 

ECZ) 

RESULT = 

UE 


GO 

TO 

40 




20 

IF 

CZ 

.LE. 

ECZ) 

RESULT == 

UE/ECZ 


IF 

CZ 

.GT. 

ECZ) 

RESULT t= 

0.0 


GO 

TO 

40 




30 

RESULT 

- 0 

.0 



40 

CONTINUE 





RETURN 

END 
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APEEKDIK C 

PROGMM LCTCSD: A XJSEP’S MAMK.I. 


Program LCYC3P calculates the nonlinear stability characteristics 
of the combustion chamber described in Pig. 3 by numerically integrating the 
system of differential equations given by Eq. (20) . Except for the term 
.-3 this equation -is the same as Eq. (12) of Ref-, -li, -whose ' 
solution is carried out by the program Lcyc3D described in detail in Appendix D 

• '• j 1 

of Ref. 11. The present computer program is Tety similar tb Program LCTC3D 
of Ref. 11 in its general structure, input and output. Hence in this user’s 
manual, only the complete listing of the present program, .along -with a precise 
description of the necessar^^ input, is given;' for details about the program 
(including input) one is referred to Appendix B of Ref. 11. 

Ho. of ' ' . . 


Cards 

Location 

Type 

Input Item 

Comments ‘ 

1 

1-5 

I 

NOUTCF 

If 0: coefficients are not 
prin|:ed out - ' 

If 1: only^ the linear coeffi- 
cients are printed out 
If 2: all the coefficients 
are printed out 


6-10 

I 

H0ZHI2 

If 0: nozzle nonlinearities 
not included 

If 1: nozzle nonlinearities 
included 

1 

1-72 

A 

TITLE 

Title used to label the 
plots 

1 

1-10 

F 

EH 

Interaction index, n 


11-20 

F 

TA.U 

Time lag, f 


21-30 

F 

H 

Time increment for numerical 
integration 


31 -Ml 

F 

TSTART 

Time at which output of 
solution begins 
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No. of 


Cards 

Location 

lirpe 

Input Item 

Comments 


41-50 

F 

TQUIT 

Time at ■which output of 
solution ends 

1 

1-5 

I 

IJTEST 

If 0 : confute transient 
■behavior 

If 1 : compute limit -cycle 

■behavior 


6-10 

I 

JMODE 

Identifies the amplitude - 
function used to test for 
limit-cycles 


11-15 

I 

UIOC 

Determines location for 'wall 
pressure maxima' and minima 





If 1 ; z; = Oj 0 = 0 ° 
If 2 : z = 0 , 0 = 45° 
If 3: z = 0, 0 = 90° 


16-20 

I 

NIEEMS 

dumber of amplitude 
functions given initial 
values 


21-25 

I 

WZ 

Determines how secondary 
insta^bility zones are 
handled 

If 0 ; all insta^bility zones 
included 

If 1: secondary zones 

eliminated 


26-30 

I 

mm 

Determines output 
If 0 : printed output 

only 

If 1 ^ EODT ^ 6 : both 

printed and plotted output; 
NOUT being the number of 
the last plot produced 


31-35 

I 

ICTYEE 

If 1 ; amplitudes selected 
to satisfy the nozzle 
boundary condition 
If 2 : amplit'udes selected 

to eliminate the exfcraneo'us 
solution 
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The next three cards are necessary only if 1 £ KOUT ^6, , 
Ho. of 


Cards 

Location 

Type 

Input Item 

Comments 

1 

1-10 

F 

YHI(I) 

Maximtim ordinate for 
pressure plots 


11-20 

F 

IHI(5) 

Maxim-Tim ordinate for 
velocity plots 


21-30 

F 

YLAB(l) 

Interval for ordinate 
labeling of pressure plots 


31"4o 

F 

YIAB(5) 

Interval for ordinate 
labeling of velocity plots 

1 

1-5 

I 

ITICI(I) 

Humber of ordinate tic 
marks for pressure plots 


6-10 

I 

ITICI(5) 

Number of ordinate tic 
marks for velocity plots 


11-15 

I 

33FIEST 

dives the number of the 
first plot produced 


16-20 

I 

HpMIT 

If 0: time -history plot 
produced 

If 1: time-history plot 

omitted 

1 

_l-5 

I 

mdplotCi) 

If 0: plot of the first 

mode amplitude not 
produced 

If 1; plot of the first 
mode amplitude is produced 


6-10 

I 

]yii)PLor(2) 

If 0: plot of the second 

mode amplitude not produced 
If 1: plot of the second 
mode amplitude is produced 


11-15 

I 

MDPL0T(3) 

If 0: plot of the third 

mode amplitude not produced 
If 1: plot of the third 

mode amplitude is produced 
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TSo. of 
Cards 


location Type Biput Item Comments 

16-20 I MDPL0T(4) If 0: plot of the pressure 

amplitude of the first 
mode not produced 
If 1: plot of the pressure 

amplitude of the first mode 
is produced 


The next card is necessary only if plot of any mode -amplitude is desired. 


1 1-10 P 

11-20 !■ 

21-25 I 

NTERMS 1-5 I 


YHB© 

XIABMD 


ITICMD 


J 


Maximum ordinate for mode- 
aaplitude plots 

Interval for ordinate 
labeling of mode-anplitu.de 
plots 

JTumber of ordinate tic 
marks for mode -amplitude 
plots 

Identifies complex amplitud« 
function 


6-15 

F 

AST 

16-25 

F 

ACT 


Amplitude of sin(u)t) terms 
in initial conditions 

Amplitude of cos (tut) terms 
in initial conditions 


The next card is necessary only if ICTIjEH - 2. 

1 1-10 F DAMP Damping factor in initial 

condition, obtained from 
linear stability analysis 
(Appendix E of Ref. Id) 

11„20 F FREQ Corresponding frequency 


in 



FORTRAN Listing 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


s 

***+jtt+***++n,+*+ pftOGJjiiM LCVC3D ++*+**♦*♦+♦*****♦♦***♦♦♦+♦**#** 

THIS. PROGH/iM CALCULATES THE NONLINEAR REHAVIOF. Of 
TRANSVEPsSE.r AXIAL.» OB COMBINED LONbI TUDINAL- TRANSVERSE 
instabilities IN A CYLINDRICAL COMBUSTION CHAMBER VI TH 
UNIFORM PROPELLANT IN<JECTION.» DISTRIBUTED COMBUSTION 
PROCESS^ AND A CONVENTIONAL NOZZLE. THE COMjEUSTION PROCESS 
IS DESCRIBED BY CROCCC'S TINE-LAG MODEL. BOTH TRANSIENT 
AND LIMIT-CYCLE SOLUTIONS ARE CALCULATED. 

THE FOLLOWING INPUTS ARE REQUIRED 

(1> THE CONTROL NUMBERS^ NOUTCF AND K0ZNL2. 

(S) THE coefficients FROM PROGRAM COEFP'SSD. 

<3) THE DATA DECK. 

NOUTCF DETERMINES PRINTOUT OF COEFFICIENTS. 

IF NOUTCF = C COEFFICIENTS ARE NOT PRINTED OUT. 

IF NOUTCF *= I LINEAR COEFFICIENTS ONLY ARE PRINTED OUT. 

IF NOUTCF = 2 ALL COEFFICIENTS ARE PRINTED OUT. 

N0ZNL2 DETEPJ<!IKE£ IF THE NOZZLE NONLI NEAFI TI ES ARE TO BE INCLUDED. 
IF N0ZNL2 c 0 NOZZLE NONLI NEARl TI LS NOT INCLUDED. 

IF NOZNL2 = I NOZZLE KONLINEARI TI ES INCLUDED. 

THE DATA DECK CONTAINS THE FOLLOWING INFORMATIONS 

TITLE OF THE RUN. 

EN IS THE INTERACTION INDEX. 

TAU IS THE TIME LAG. 

HIS THE INTEGRATION STEP SIZE. 

TSTART IS THE TIME AT WHICH OUTPUT STARTS. 

TOUIT IS THE TIME AT WHICH COMPUTATIONS ARE TERMINATED* 

NTEST IS TASK CONTROL NlKBEFi 

IF NTEST = O COMPUTE TRANSIENT BEHAVIOR. 

IF NTEST = i COMPUTE THE LIMIT-CYCLE BEHAVIOR. 

JMODF IS THE MODE- AMPLITUDE USED TO TEST FOR LIMIT-CYCLES. 

NLCC DETERMINES THE LOCATION OF THE WALL PRESSURE MAXIMA' 

AND MINIMA; 

IF WLOC c 1 LOCATION I S Z » o, THETA = 0 DEGREES. 

IF NLOC = 2 LOCATION IS Z = 0* THETA = A5 DEGREES. 

IF NLOC = 3 LOCATION IS Z = 0* THETA a 90 DEGREES. 

NTERMS IS THE NUMBER OF TEIJ<?S GIVEN INITIAL VALUES. 

NPZ DETERMINES ROW SECONDARY STABILITY ZONES (PHANTOM- 
ZONES) ARE HANDLED. 

IF NPZ - 0 PHANTOM ZONFS ARE RETAINED. 

IF KFZ = 1 PHANTOM ZONES ARE ELIMINATED. 

NOUT IS THE OUTPUT CONTROL NUMBER. 

, IF NOUT = 0 printed OUTPUT ONLY. 

IF NOUT > 0 BOTH PRINTED AND PLOTTED OUTPUT* NOUT 
DETEi^INES THE KUMEEF. OF THE LAST PLOT 
PRODUCED. 

ICTVPE IS THE INITIAL CONDITION CONTROL NlMEERt 
. IF ICTYPE = 1 AKFLITUEES SELECTED TO SATISFY 
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THE KOZZLE BOUNtAHY CONLITION. 

IF ICTYFE = 2 AlOFLITBDES SELECTEE TO ELIMINATE THE 

EXTRANEOUS SOLUTION. 

DATA FOB SETTING UP PLOTS % 

YHI(l) IS THE MAiilNW^ ORDINATE FOE FBESSUBE PLOTS* 

YHX(5> IS THE MAXIMUM ORDINATE EOH VELOCITY PLOTS* 

NOTEt THE ORDINATE SCALES FOB PRESSURE AND VELOCITY PLOTS 
ABE SYMMETRIC ABOUT ZERO* 

YLAB IS THE INTERVAL FOR ORDINATE LABELING FOR ABOVE PLOTS. 

ITICY IS THE NUMBER OF ORDINATE TIC MARKS FOR ABOVE PLOTS. 

NOTE; ITICY SHOULD EE NEGATIVE FOB PRESSURE AND VELOCITY PLOTS 
TO OBTAIN CENTERLINE. 

NFIRST IS THE NUMBER OF THE FIRST PLOT PRODUCED. 

NOKIT DETERMINES VHETHER AMPLITUDE PLOT IS PRODUCED 
IF NOMIT = 0 AMPLITUDE PLOT IS PRODUCED. 

IF NOMIT 1 AMPLITUDE PLOT IS OMITTED. 

MDFLOT DETERMINES IF THE PLOT OF THE MODE-AMFLI TUDE IS RE&UIHED* 
IF KDPLOT “ 0 PLOT NOT REQUIRED* 

IF MDPLOT = I PLOT REQUIRED* 

YHIMD IS THE MAXIKLSO ORDINATE FOB AMPLITUDE PLOTS* 

YLAEMD IS THE INTERVAL FOR ORDINATE LABELING OF AMPLITUDE PLOTS* 
ITICMC IS THE NUMBER OF ORDINATE TIC MAl2iS. 

NOTE: ITICMD SHOULD EE NEGATIVE TO OBTAIN THE CENTERLINE* 

INITIAL AMPLITUDES OF F-FUNCTIONS CREMAINING CARDS) 

AS<J> IS THE AMPLITUDE OF THE SINE TERM. 

ACCJ> IS THE fWPLITUDE OF THE COSINE TEPK* 

DAMP AND FREO ARE THE DAMPING COEFFICIENT AND THE FREQUENCY FROM 
THE LINEAR STABILITY PROGRAM* 


PARAMETER 

COMPLEX 

COMPLEX 

DIMENSION 

1 

2 

3 

4 

5 

6 
7 
B 
9 
1 
2 
3 


~KX=5# KX2elO* MXAeSOi MXSSQelOO 

YNOZ(MX)> B<MX)* Cl* Cg* C3* CPHITCMX)* CSW* A 
GNOZ(MX)* CAXl* Cl 

L<NX>* N(MX)* S(MX)* NAME<MX)> AS<MX2)* AC<MX25* 
U<250>MX4)* Y(MX4)/ FZt4*KX4>* YP<MX4)» UZ<KX4)* 
CP<4#MX2.MX2)* FRGKMX2)* IMPHMX23* UMAXCS005* 

ZC6). ANGLEC6)* IHETA(6>* CFT<6*MX8)* YKMX2)* 
CFTHC6/MX2)* CFZ<6*MX2>* PRESS<6)* AXVELO)* YB<MX2)* 
TPLOT<500)> YPL0TC6* 500)* DWMYT<500)> DUMMYY< 500>* 
IBUFOOOO)* ITTC 4 )* ITYK 7 )* ITY 2<73* nY3<7)ji 
1TY4(7 )* ITY 5<6 )* TAUCUT<MX2)* 1TY6<8>* UAVGUOO)* 
1TP<3). TITLE(12). FRSC500)* TK500)* FMAX(500>* 
TIMAX<500). YL0C6)* YHK6). YLAB<G)* ITlCYCfe)* 
KFREOCKX)* VKPCMX)* AAt4)* UPLOTCMX* 500)* FRlTtSOO)-* 
MDPL0T(4)* MTITLK4)* MTITL8(4)> MTITL3 CA)* 

*"^1TL<4 )* PRTITLC5) 





c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


COMMON 

1 

2 

COMMON 

COMMON 

COMMON 


RVCMXa>4)» CC4*MX2>MX23» D<MX2>KX2£Q) > 
KPMAX<4#KX2)> iC(4>MX2^MX2)^ KpGMAX<MX23i 
IDF<WX2>MX2S03 j I D0tKX2*MX2S&) 

/BLKg/ KCKX3J. NS(MXJ> £J(KX>> B 

/BLK3/ NOMAXj KLMAX> 6AMKA» COEFOjMX23 

/NLTEFM/ N0ZNL2rf EXThA(MX£jA) 


DATA 

1 

2 

3 

A 

5 

6 
7 
5 
9 
1 
2 


I TT/*’ DIMENSIONLESS TIME# T"/, 

I TYl/” INJECTOR FEESSORE PERTURBATION# THETA = O’*/# 
nY2/**IKJECTOR PRESSURE PERTURBATION# THETA * AS"/# 
ITY 3/ "INJECTOR PRESSURE PERTURBATION# THETA s= 90"/# 
ITYA/"NOZZLE PRESSURE PERTURBATION# THETA c 0"/* 
ITYS/”NOZZLE AXIAL VELOCI TV# THETA = 0*'/# 

I TY6/**K0ZZLE B»C* ( RE<-GAHMA*Y*PHI T) > AT THETA = 0" /# 
1 TP/"PRE£SURE PEAKS**/ 

MTITLl/" AMPLITUDE OF IT MODE*'/ 

KTITL2/"AMFLITUDE of 2T MODE"/ 

MTITL 3/ "AMPLITUDE OF IR MODE**/ 

PHTI TL/*'PRE£SURE AMPLITUDE OF IT MODE"/ 


LAST «= 250 
ERR »= 0«001 
TDEL *: 10*0 
NPT = 0 
AACn = 0*0 
AA<2J e 0*5 
AA<3> c 0*5 
AA<A> = I#0 
PI = 3-1415927 

READ C5#5003> NOUTCF# N02KL2 


♦♦**^=%**+**^* COEFFICIENT INPUT SECTION ******++:^.,c*++>^+ + H=******=(c«* 

THIS VERSION OF LCYC3D READS THE COEFFICIENT DATA FROM 
A FASTPAND FILE GENERATED BY PROGRAM C0EFBS3D- TO BEAD 
THIS DATA FROM CARDS# USE READ <5#XXXX) INSTEAD OF 
READ C9#XXXX) IN THIS SECTION. 


INPUT OF MOTOR PARAMETERS AND NIMBER OF TERMS- 

READ (9# 5001) GAMMA# UE# ZE# 2C0ME# NDROFS# NJKAX# NOZNLl 

VRITF (6#600I> GAMMA# UE# ZE# ZGOKB# KJMAX 

IF CK DROPS -to. OJ UHIIE C 6# 6030) 

IF CNDHOPS -EO. 1) VRl TE C 6# 6031) 

IF <N0ZNL2 .EO. 0) Iv’RI TE C 6# 6032) 

IF <K0ZNL8 -EO. 1) VRI TE C 6# 6033) 

NU = 2 ^ NJMAX 
JMX = NJMAX/-2 
RLD s= 0-5 * ZE 


WRITE (6# 6002) 

C 

C INPUT OF DESCRIPTION OF SERIES EXPANSION. 

DO 10 K = 1# JMX 

READ <9#5002> NJ# LfNJ)# M<NJ># K<NJ)# NSCNJ)# SCNJ)# SJCNJ)# 

1 ki Aim V • 1 % 


ll4 
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WRITE C6>6003) NAME(NJ)j NJj LCNJ)« MCNJJj N(NJ>« 

1 SCNJ}> SJCMJ} 

10 CONTINUE 

WRITE C6>6010> 

DO 15 K e I, JMX 

READ C9>5010> J* YN0ZCd)> B(J) 

WRITE C6>6015> J» YN0Z<d># ECJ) 

NJ e C2 * JJ - 1 
YRCNJJ »= REALCYN02(d>> 

YICNJ> *= AIMAGCYNOZtU>> 

YRCNJ+n = YRCNJ) 

YICNJ+1) YICNJ) 

15 CONTINUE 

IF (NOZNLl .NE. 1>,G0 TO 815 
WRITE <6#603A> 

DO 820 K e 1, dWX 

BEAD C9#50in J* GNOZCdJ 

WRITE <6*6035> GNOZCd) 

820 CONTINUE 
815 CONTINUE 

CALCULATE THE NUMBER OF TYPES OF LINEAR COEFFICIENTS. 

NCOEFF «=^ A 

IF <N0ZKL1 .EQ. 15 NCOEFF = 5 
NCFWl e NCOEFF -I 

ZERO LINEAR COEFFICIENT ARRAYS- 
DO 20 KC t: NCFMI 

DO 20 Nd 1* MX2 

DO 20 NF t 1. MX2 

C<KC>Nd.NP5 c 0.0 
CP<KC.Nd.NP5 «= 0-0 
20 CONTINUE 

ZERO NONLINEAR COEFFICIENT ARRAY - 
DO 30 Nd « 1* MX2 
DO 30 NPQ *=_l,_MX2S0 
DCNd^NPO) t= 0-0 
30 CONTINUE 

INPUT OF LINEAR COEFFICI ENTS- 
DO AO KC = KCFMl 

READ C9> 50035 KMAX 

IF CNOUTCF -GT- .05 WRITE <6. 600A) KC. KMAX 
IF IKMAX -EQ* 05 GO TO AO 
DO A5 K = 1> KMAX 

READ (9.500A5 Nd> NPj CP<KCjNd*NP) 

IF (NOUTCF *GT. 05 WRITE <6^60055 KC> Nd» Np> CP<KC.Md^NF5 
AS CONTINUE 
AO CONTINUE 


INPUT OF N0NLI^3EAR COEf FI Cl ENTS. 
READ <9# 50035 NLMAX 
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li- tNOUTCF .E0. 2) WHITE <6*6006} NLMAX 
IF <NLbi?AX .EQ. 0} GO TO 50 
CO 52 NJ = 1* MX2 
KPQMAX<Nd) c*'o 
52 CONTINUE 

DO 55 K *= I, NLMAX 
BEAD C9*5005) NJ* NF* NO* DT 

IF (NOUTCF .EG* 2} WRITE <6*6007> NJ# NT* NQ* DT 
KFGNAXCNd) »= KPGKAXINd) + 1 
KPO *= KFGMAXINJ} 

IDF<NJ,KPQ5 = NF 
1D0<NJ*KP0) = NO 
D(NJ*KPQ) = DT 
55 CONTINUE 

50 CONTINUE 

♦ *****.^++**+* PRESSURE COEFFICIENT SECTION »♦*»♦**#♦♦♦♦**#++♦♦*** 

CALCULATE SPATIAL COORDINATES FOB PRESSURE COMPUTATION* 

DO 51 NFRES = 1* 3 
Z(NPRES) <=0*0 
RTHETA = NPRE'S - 1 
ANGLEINPRES) *= RTHFTA * 45*0 
THETACNPRES} = RTHETA * FI/4.0 
ZCNPBES + 33 ■» ZE 
ANGLE< NFRES + 3) = ANGLE<NPRES3 
THETA(NPFiE£ + 3> = THETACNPKES3 

51 CONTINUE 

CALCULATE COEFFICIENTS FOR PRESSURE TIME HISTORIES. 

DO 53 NPRES = 1* 6 
DO 53 J = 1, dKX 
NP = (S * dJ - I 
Z1 c Z<NPRES} 

ANG = THETACNPRES) 

CALL FHICFSCd*Zl*PNG*Cl*C2*C3) 

IF CNFRES .E0. 4) CPHITCd} « Cl 
CFT<NPRES.NF> t= REALCCD 
CFTCNPBES*NP+ n = -AIMAGfCl) 

CFTH<NFRESi NP) = REALCC2) 

CFTH<NPPE£*NP'f- 1) b “AIMAG(C23 
CFZCKPRES*NF> REAL<C3> 

CEZCNPRES*NF+1> = -AI MAG CCS 3 
53 CONTINUE 
C 

Cl = co.c*i»o> 

CAXl e GAMMA * CCOSHCCI + B<1) * ZE) 

CAXIR e RFAL<CAXI3 
CAXII = AIMAGCCAXI) 

C 

C OUTPUT OF' COEFFICIENTS FOR PREISSUBE TIME HISTORIES. 

WRITE <6*6020) 

DO 56- NPRES «= 1* 6 
• WRITE <6*6014) 

DO 56 d c I* NdMAX 
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VFilTE <€j6081> Z(NFBES)> ANCLECNFEES)^ 

1 CFKNPRES* J)> CFTH<NPBES*0)^ CFZ < NFRES^ J) 

56 CONTINUE 
C 

C *♦»♦♦♦♦♦♦♦♦♦♦ data input section ♦^♦♦♦♦♦♦♦♦************+***=<'*»***’^ 

c 

READ C5>5000) TITLE 
C 

C ZERO INITIAL VALUE AND FREQUENCY ARRAYS. . 

5 DO 57 K K 1# NJKAX 
AS(K> ?= 0.0 
ACCK) *=0.0 
FROKK> *: 0*0 

57 CONTINUE 
C 

C 

C READ COMBUSTION AND CONTROL PARAMETERS. 

READ (5>5006jp END = 300> EN*. TAU# B* TSTART* TGUIT 
C 

C READ CONTROL NUMBERS. 

READ <5»5008J NTEST# JMODE» NLOC# NTERMS# NFZ» NOUT# ICTYFE 

JMODE «= ce ♦ JM0DE3 - 1 

JFMODE = OMODE + NJKAX 

IF <N0ZNL2 .NE.’lJ GO TO 825 

FEES *= SC1> 

KFREOC n *= I 
KFRE6<2) = 2 
KFREOO) ** 2 
DO 830 K e 1 j JKX 
VKPCJ) *: FKEG * KFEEQtJ) 

830 CONTINUE 
825 CONTINUE 
C 

IF CNOUT .GT. 05 NFT & 1 
IF (NOUT .EQ» 0) GO TO 9 
C BEAD DATA FOR SETTING UP PLOTS* 

READ <5*50095 YHIC15* YHK55* YLABC15* YLAE<55 
READ <5*50085 ITlCY(l)* ITICY<55* NFIBST* KOKit 
READ <5# 50145 MDPLOT 
MDFLTL =* 0 
DO 3S0 K c I* JMX 
MDFLTL = MDFLTL + MDFL0T<K5 
320 CONTINUE 

IF (MDFLTL *EQ. 0> GO TO 9 
READ <5*50155 YHIMD* YLAEMD* ITICMD 
YLOMD - - YHIMD 

♦*♦♦♦*+♦**** INITIAL AMPLITUDES SECTION ♦♦♦♦**♦**♦♦♦♦*♦♦♦+*♦*♦**** 
9 DO 58 K = I* NTEBMS 

p 

INPUT INITIAL AMPLITUDES FOR F-FUNCTiONS. 

READ (5*50075 J* AST* ACT 
NJ = <2 * J> “ 1 
AS(NJ5 c AST 
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C 

C CALCULATE FREQUENCY AND DAMPING. 

IF <ICTYPE .EG. 23 GO ' TO 58A 

RL «= LC03 

AX t= RL * PI/ZE 

AX SO = AX ♦ AX 

SSQ SCJD * SCJ3 

FR0KNJ3 c SQRTCSSe + AXSQ3 ■ 

ENFKNJ) = 0.0 
GO TO 586 
584 LONG = LCJ) 

SWN = SCU) 

READ (5> 5099) DAKF*FREQ 

DNP1<NJ> = daf:f 

FRCKNJ3 = FBFO 
586 CONTINUE 

FEQICNJ+13 = FRQ1CNJ3 
DMPKNJ+13 c DMF1CNJ3 
C 

IF flCTYPE .EG. 2) GO TO 582 
C CALCULATE INITIAL AMPLITUDES FOE G-FUNCTIONS* 

C 

IF CFROKNJ33 58# 58# 58 1 
581 GYRU - GAKMA*YR(Nd3*UE 

GYIF = GAMKA^YI CNd>*FEGl(NJ) 

GYEF = 6AKMA+yR(NJ>=^FRGlCNJ> 

GYIU = 6AMMA*YI (NJ3*UE 
C 

NFBES c 4 

IF tNS<d) -EQ. 1) NFRES = 6 
C 

A1 = Cl.O + GYRU)*CFZCNFRES#KJ+ 1) 

I - 6YI F+CFT(NPRES#NJ+ 1 ) 

AS s= GYRF*CFTCNPhES#NU+ 1) + GYI U*CFZ<NpfiES# KJ+ 1 3 
A3 = -<1.0 + GYRU3 + CF2CNPEES>Nd3 +. GYIF*CFTCNFRES#NJ3 
A4 = GYRF«CFTCNFRE£#NJ3 + GYI U4=CFZCNPRES# Nd) 

C 

DET - A1*A1 + A2*A2 
IP.CDET .LT. 0.00000013 60 TO 583 
III *= A3+ACCNd3 - A4+ASCKdD 
R2 = -A4+AC«NJ) - A3+ASCN03 
C 

ACCNJ+13 = CRl+Al + B2*A23/DET 
A5CKJ+13 = -CR2+A1 - R1*A23/DET 
GO TO 58 

583 ACCNJ+13 = -ASCNd) 

- ASCNJ+13 a ACCNJ3 
GO TO 58 
C 

582 AEG- = FROlCNd) ^ TAU 
FSIN = SINCARG3 
FCOS a 1. - C0SCARG3 
FSO = FRQICNJ) * FRCHNJ) 

DSC a DMPKNd3 * LKFlCNd> 
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Al e DSQ - FSQ + EMFICNJ) * < CP C N J* M J ) 

1 - EN * CF(3>NJ*NJ) * FCGS) 

2 + Eh.* CF<3>NJiNJ> * FRQICNJ) * FSIN 

3 + CPC1,NM,NJ) 

AS = CS.O * I»?FICNJ> + CPCSjNJ^NJJ 

1 - EN * CPO/NJ/hJ) * FCOS) * FRGICNJ) 

2 - EN * CFO..NJ,KJ> * EMPKNJ> + FSIN 

A3 «= CPtSjNJjNd+ 1> * EMPICNJ5 * CPCl>Nd*KO+l> 
AA = CF<2*Nd»NJ+l > * FRS1(NJ> 

DEN *: A3*A3 + A/|*A4 

IF (DFN .LT. 0.0000001> GO TO 585 

R1 = A1 + A3 +AS*A/( 

R2 s Ai*A4 - A2+A3 

AC<Nd+l) = <-hl*ACCNd) + R8+AS(K0> )/DEN 
ASCNJ+1) = -<RS*AC(NJ> + KJ *ASCNO) > /PEN 
GO TO 58 

585 AC(NJ+1> s -AS(NJ) 

AS(NJ+1> e AC(NJ> 

C 

58 CONTINUE 


OUTPUT OF INITIAL AMPLITUDES. 

WRITE (6>6016) 

DO 59 0 J = 1. WJKAX 
IF CA5CJ3) 59 1 1 59 2. 591 
592 IF (ACtJ)) 59 1. 590. 591 

591 WRITE <6.6017) J. DMPICJ). FP.OICJ). ACCJ). ASCJ) 

590 CONTINUE 

IF CNTEST .EG. 0) WRITE (6.6025) 

IF (NTEST .EO. 1) WRITE (6.6026) 

IF (NPZ .EO. 1) WRITE (6.6028) 

IF (NOUT .GE. 1) WRITE (6.6027) 

*4:*****«4c«**« LINEAR COEFFICIENTS SECTION ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦H:**** 

DO 59 KC » 1. NCFMl 
DO 59 Nd - 1. EX2 
KPMAX(KC.Nd) t= 0 
59 CONTINUE 

IF (NPZ .EO. 0) GO TO 605 
DO 602 d = 1. dMX 
KJ = (2 * d) - 1 
RL = L(J) 

AX = RL * FI/ZE 
AXSO c AX * AX 
SSC = S(d) * S(d) 

OMEGA = SCBT(SSG + AXSQ) 

TAUCUT(Nd) = 2.0 * FI /OMEGA 
TAUCUT(Nd+l) t= TAUCUT(Nd) 

602 CONTINUE 

DO 604 NJ = 1. NdKAX 
DO 604 NP = I. NdKAX 
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IF <T/iU *GT. TAUCUTCNFJJ CFC3 #NJjNF> = 

604 CONTI WUE 

COMPUTE LINEAR COEFFICIENTS FOR GIVEN VALUES OF EN AND TAU* 

605 DO 60 NJ *: NUMAX 

DO 60 NF = 1» NJKAX 
CT = CP<l#NJ/NP3 

IF <CT) 61* 68* 61 
61 KPMAX<1*NJ) = KPMAXC1*NJ3 + 1 
KP = KFMAXC1*NJ3 
IC<1*KJ*KP) = NF 
CC1*NU*KP) = CT 

68 CT = CPf8*N'J*NP) - EN+CFC 2*N0*KP) 

IF CCT> 63* 64* 63 

63 KPMAX(8*KJ) *= KFMAX(a*WJ5 + 1 
KP = KPMAX(2*NJ3 
ICCS*NJ*KP) c NF 
C<8*NJ*KP) s= CT 

64 CT = EN ♦ CPC3*NJ*NP> 

IF CCT) 65* 66* 65 

65 KPKAXC3*KJ) = KFMAX(3*NJ> + 1 
KF KPMAX(3*NJ) 

ICC3*NU*KF) NF 
Cf3*NJ*KP) = CT 

66 IF CN0ZNL8 *NE. 1) GO TO 60 
CT = CP<4*Nd*NF3 

IF CCTi 67*60*67 

67 KPMAX<4*NJ) KPMAX<4*Nd> + 1 
KP c KPWAX<4*NJ) 

ICC4*NJ*KP) e NP 
C<4*Nd*KP) = CT 

60 CONTINUE 

STEP- SIZE COMPUTATION ***^***:¥****************^:fi**yt^=t^ 

NDIV = 1*0 + TAU/H 
RN = NDIV 
H = TAU/FN 
R6 e H/6*0 

4: *4; INITIAL VALUES SECTION ******4:* *+***++ + ** *****+*+ + + + 

VHITE C 6* 6008) EN* TAU* GAMMA* UE* RLD 
VRITF C 6* 6009) 

V/RITE C6*602S) CANGLE(J)* J s= 1*6)* (ANGLECJ)* d = 1,3> 

WRITE (6*6012) 

NFl = NDIV +1 ■ 

DO 70 I 1* NPl 

NSTEP t= 1 - NPl 

RSTEP c NSTEP 

TIKE e RSTEP * H 

TKI) = TIME 

DO 75 d = 1* NdKAX 

dP = d + NdMAX 

IF <AC(d)> 7S1* 753* 751 
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753 IF <AS<J)5 75J* 758* 751 
752 IKI*J5 e 0*0 
U<1*JP> B 0*0 
GO TO 75 

751 ARG *= FRQKJJ ♦ TIME 
ESIN <= SI N< ARC 7 
FCOS *= COSCAKGJ 
FEXP «= EXP<rMPlCd7^TIME) 

U<1*0) = <AS<J>*FSIN + AC<0)4FC0S5 + FEXF 

UCI*JP) e <<A£<J) * FCOS) - <AC<J) ♦ FSIN)) * FR8HJ) # HXF 
1 + DMPlCd) * UCI*J) 

75 CONTINUE 

C CALCULATE INITIAL VALUES OF PRESSURE AND VELOCITY. 

DO 704 NPRES *= 1* 6 
DO 708 U = I* NJMAX 
C0EF<1*J) *= CFT(NPEES*J) 

C0EF(2*d> = CFTH( NPRES* U) 

C0EF<3*0) « CFZCNPRES*<J) 

702 CONTINUE 

DO 703 J » I* NU 
Y<J) « U<I*0) 

703 CONTINUE 
UBAR =0.0 

IF CNPRES -GT. 3) UBAR *= ' UE 
= 0*0 

IF < (NDROPS-EO. U .AND. (NFRES.LT. 4) ) UKS = UE/(ZE+ZC0ME> 

CALL FRSVEL<UBAR*UMS*Y*P*VTHi UZ) 

FRESS< NPRES) = P 

IF' (NPRES .GT. 3) AXUEL(NFRFS - 3) = VZ 

704 CONTINUE 

PRS(I) = FBESSINLOC) 

C 

C CALCULATE INITIAL VALUES OF NOZZLE B.C. 

CSUM = (0. 0*0.0) 

DO 710 J = 1* JMX 

JF = NUMAX + (2 ♦ J) - 1 

FT = Y(JP) 

GT = y<JP+l) 

A = CKFLX(FT*GT) 

CSUM = CSUM + VNOZCJ) * CFHITCJ) + A 
710 CONTINUE 

SUM = REAL C CSUM) 

YFHI = -GAMMA * SIM 

URITE C6*601l) NSTEF* TIME* (PRESSCJ)* d = 1*6)* 

I (AXVELCd)* 0 = 1*3)* YFHI 

70 CONTINUE 
C 

WRITE <6*6008) EN* TAU* GAMMA* UE* RLD 

WRITE (6*6088) CANGLECd)* J = 1*6)* (ANCLE(J)* J = 1*3) 

C 

c INITIALIZE CONTROL NUMBERS ****** +++ + *+*>l<********++* 

C 

LINE = 8 

K = 0 • 

MAXNO = 0 



o o non 


KAXP = 0 

IF fNOlJT •EG* 0> GO TO 100 

OFLOT = 0 

TWIN = TSTi!4FT 

TMAX = TSTAET + TDEL 

YL0< n = -YHKlJ 

DO 90 J Sf M 

Yhi(j) = YHicn 

YLOC J) 5= YLO< 1) 

VLAB<J) = YLAB< 1) 
niCY(J> t= IT1CY<1) 

90 CONTINUE 

YLOtS) c “YHIC5> 

YHKe) c YhICS) 

VL0<6) = YL0C5) 

YLAE<6) = YLAB(5) 

ITICY<6) = ITICY<5> 

♦*****♦**♦*♦♦ NW^EBICAL CALCULATIONS SECTION ++*****+*+**■**++** 
100 I = NPl 


RiaJGE-KUTTA INTEGRATION SCHEME. 

105 NSTEP c Cl - NFl + (LAST - NF 1 ) * K> 
KSTEP = NSTEP 
TINE = PSTLP * H 
TICI) = time 
DO 110 O = Ij NJMAX 
OP J + NOMAX 
RVCO^IJ c UCI-NDIV^OF) 

RVCO, 4) = UCI-NDIV+IjOP) 

KUCJ^SJ = 0.375+RVC0j 1) + 0»75*RVCJ>4) 
RVCJ, 3) = RVCJ, 2) 
no CONTINUE 

IF (NOZNLE .NE. 15 GO TO 835 

DO 840 II = 1,4 

TZ = TIME + AACII5+H 

DO 840 d »= 1,JMX 

00 DD = S*J - 1 

OEVEN c 2«=J 

EXTEACdODD,II ) = COSCVKF (0)*TZ) 

EXTRA < OEVEN, I I ) = SlNCWKFCOl^TZ) 

840 CONTINUE 
835 CONTINUE 

DO 120 0 = 1, NO 
Y<d5 = UCI,d5 
120 CONTINUE 

CALL RHSCNU, 1,Y,YP> 

DO 130 d =■ 1, NU 
FZCl,d) c YPCdl 
130 continue 

DO 140 XI = 2,4 
DO 144 d = 1, NU 

UZCd5 = YCJ) + AACII5 + B * FZClI-I,d) 
144 CONTINUE 


0. 125*UCI-NDI V+2,dP5 
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c 

c 


C^lLL RHSCNU*11>UZ*YP) 

DO \AZ 0 *= NU 
FZai^O) *= YF<J) 

148 CONTINUE 
140 COKTINUE 

DO 150 J «= 1* NU 

U<I+l>J) t: TCJ) + (FZ< 1j J)+2.0*tFZ<2*U>+FZ«3^0) ) + 
150 CONTINUE 


FZ<4»J>) * H6 


CALCULATE PRESSURE TIME HISTORIES* 

DO 154 NFRES >=1/6 
DO 15S J = 1/ NdMAX 
C0EF(1>J) *= CFT<NPRES>J) 

C0EF(2/d) = CFTHCNPRES* J> 

C0EFC3/J> *s CFZCNPRES/J5 
152 CONTINUE 
UBAR *= 0*0 

IF CNPRES *GT. 3) UEAR c UE 
UMS = 0*0 

IF CCNDRGPS.EO. 1 J *AND* CNPRES-LT* A) ) UMS «= UE/CZE»ZCOME> 
CALL PRSUEL< UEAR/ IW S# Y/ P/ UTH/ VZ ) 
presscnfp.es > *= P 

IF CNPRES *GT* 33 AXVELCNPBES - 33 = UZ 
154 CONTINUE 

PRSCI3 = FRESSCNL0C3 


C 

C CALCULATE VALUES OF NOZZLE B*C* 

CSUM = <0*0/0*03 

DO 650 J = 1/ JMX 

dP = NdMAX + <£ ♦ d3 - I 

FT = Y<dP3 

6T e YCdP+13 

A CMPLXCFT/6T) 

CSUK = CSUM + YNOZCd) * CPHITCd> * A 
650 CONTINUE 

SW = EEALCCSUM) 

YPHl *= -GAMMA * SUM 


C 

C 

C 

C 


C 

c 

c 


DETERMINE MAXIMA AND MINIMiA OF PRINCIPAL MOPE” AMFLI TUDE 
FUNCTION FOR USE IK DETERMINING LIMIT-CYCLE BEHAVIOR* 

IF <UCl/dFM0DE3 UC 1+ 1/ dPMODE) 3 170/ 170/ 160 

170 PDEN = UCI/dPMODE) - Ut 1+ 1/ dPMOCE) 

IF (PDEN3 171/ 160/ 171 

171 FP = UCI/ JPMODE)/PDEN 


PA = CPP - 1*03 * PP ♦ 0*5 
PB= 1*0- <FP* PP> 

PC = CPP + 1*03 "tf FP * 0*6 
MAXNO = MAXNO + 1 

UMAX(MAXNO> = pA*UC I- 1/ JM0DE3 + PB*UC 1/ JM0DE3 + 
IF CMAXNO *GE. 500) 60 TO 850 

160 CONTINUE - 


PC*U(I+ 1/ JM0EE3 


DETEFiMINE MAXIMUM AND MINIMUM PRESSURE AT LOCATION SPECIFIED 
BY NLOC. 
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DPL c FRS(I> - PRSCI-1^> 

DPS = PhS(I-n - FRSCl- 2 ) 

IF CDPL+DPS) 113 -, 173j 175 

173 PMUF' = FhSCI-e) “ PRSCl) 

PDEN = 2.0 + <FPS<1-2J + PRSCl > - 2« 0#PRS< I - 1 > ) 

IF CPDEWJ 174^ 175j 174 

174 PP = FNUM/PDLK 

FA e CPF - 1.03 * PP ♦ 0.5 
PB s 1.0 - CFF ♦ PF3 
PC = CPP + 1.03 * FF * 0.5 
KAXP = MAXP + 1 

PWAXCMAXPl = FA+PRSCI-23 + PB+PBSCI-13 + PC* PRSCl 3 
TIMAXCMAXF3 t= TKI-l) + FP*R 
IF CMAXP .GE. 5003 GO TO 250 

175 CONTINUE 

IF CNTEST .EQ. 13 60 TO 155 

IF CTIME .LT. TSTAHT) 60 TO 155 

IF CCNOUT .E0» 03 .OR. CKOUT .GT* 633 60 TO 156 

T-IKE HISTORY PLOTTING SECTION *+:i<***^*+4=++**=K**=*=*>i:** 
IF (TNAX .6T. TQUIT3 GO TO 156 

IF ( C TIME .GT. TKAX3 .OR* COPLOT »GE« 50033 GO TO 1000 


OFLOT = OFLOT + 1 
C 

C FILL TIME ARRAY FOR PLOTTING. 

'TFL0TCJPL0T3 = TIME 
C 

C FILL INJECTOR PRESSURE ARRAYS FOR PLOTTING C THETA = 0^ 45> 9 03 

DO 1001 J = 1>3 
YFLOTC J.OPLOT) = FRESSCO) 

1001 CONTINUE 
C 

C FILL NOZZLE PRESSURE ARRAY FOR PLOTTING C THETA = 0) 

YFLOTCAi JFL0T3 = PBESSC43 
C 

C FILL NOZZLE AXIAL VELOCITY ARRAY FOR PLOTTING C THETA = 03 
YFL0T(5^JPL0T3 = AXVELCl) 03 

C 1 

C FILL NOZZLE B.C. ARRAY EOP PLOTTING CTHETA =03. 

YPLOTCG^ JPL0T3 = YPHI 


C 

C 


32S 


IF CMDPLTL .EQ. 03 GO TO 156 


FILL MODE AMPLITUDE 
DO 322 J = 1^ JKX 
IF CMDPL0TCJ3 
J12 = 2+J - 1 
UPLOTCd^ JPLOT) 

CONTI NUE 


ARRAYS FOP. PLOTTING. 


-EQ. 03 GO TO 322 
= UCI>J123 


JlTl = 
J1T2 = 


NJMAX 

NJMAX 


1 

9 


page ts 
OP POOR QUAurr 


l2h 



oo oo r> o oo no r> no no o ooo 


PRITCJFLOT) c C^^XIR*UC1^ Jim - CAXII + UCIj J1TS5 
C 

GO TO' 156 
C 

1000 NUM *= JFLOT 

PLOT TIME HISTORIES. 

DO 1020 MPLOT <= NFIFST^ NOUT 
JFLOT = 0 

ASSIGN PLOTTING PARAMETERS. 

YMIN = YL0<NFL0T5 
YMAX »= YHICNFL0T5 
NTICY = ITICYCNFL0T5 
EELY’ = YLAECNFL0T5 

ELIMINATE POINTS THAT ARE OUT OF THE ORDINATE RANGE. 

10 1010 J = 1# NOW 

IF CCYFLOTINFLOT# JJ *LT. YKIN) -OR- < YFL0T<NFL0T^ J) .GT« YMAX)) 
1 GO TO 1010 
JFLOT JFLOT + 1 
DUMMYT< JFLOT) <= TPLOT<J) 

DUMMYYC JFLOT) = YFLOT<NPLOT> J) 

1010 CONTINUE 

IF CJPLOT .EQ. 0) GO TO 1020 
GO TO t 1011# 1012^ 1013* 1014J 101 5^ 1016)/ KFLOT 

PLOT INJECTOR PRESSURE AT THETA = 0 DEGREES. 

1011 CALL GRAPHStlEUF^SCOO^^jOFLOTiSUNTlCY^IKAX-fYKAX/TOIN^YMN# 

1 1TT>ITV1#21# A1#DWMYT#DUHMYY#2.0#DELY#T1TLE) 

GO TO 1020 

PLOT INJECTOR PRESSURE AT THETA c 45 DEGREES. 

1012 IF (MCJKODE) .£©. 0) GO TO 1020 
CALL CRAFHSdEUF# 3000# 4# JFLOT# 51#NTI CY# TM AX# YMAX# TMIN#YM1N# 

1 ITT#1TY2# 21# 42# LUMMYT# DUMMYY# 2. 0# DELY# TI TLE) 

60 TO 1020 

PLOT INJECTOR PRESSURE AT THETA « 90 DEGREES. 

1013 IE <M(JMODE) -EQ. 0) 60 TO 1020 

CALL GBAFHSCIEUE# 3000# 4# JPLOT# 5 l#KTl CY# TMAX#YMAX# TMIN#YMIN# 

1 1TT#ITY3#21#42# DWMY T# DUKMYY#'2. 0# DELY # TI TL E ) 

GO TO 1020 

PLOT NOZZLE PRESSURE AT THETA = 0 DEGREES. 
lOlA CALL GRAPKS<IBUF#3000#4# JFLOT# 51#NTICY>TKAX# YMAX# TWIN# YMIN# 

1 ‘ ITT#ITY4#21#39#DLWMYT#DUMMYY#2.0#DELY# TITLE) 

GO TO 1020 

PLOT NOZZLE AXIAL VELOCITY AT THETA «= 0 DEGREES. 

1015 CALL GRAFHSUBUF# 3000# 4# JFLOT# 51# NTICY# TMAX#YMAX# TMIN#YMIN# 

I I TT# ITY 5# 21# 32# DUMMY T# DUMMYY# 2.0# DELY # TI TLE) 
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GO TO 1020 


C 

C PLOT KOZZLE B.C. OT THETA 0 t-EGP.EES# 

lOlG CALL GKAPHS(IBUE'>3000» Aj«JFL0T> SIjNTIGY^ TMAX» T MAK» TMN/YMIN^ 

I I TT> 1 TY6j 2 OUMKYT# ClJKMYY^ 2* 0* EELY> TI TLE> 

C 

1020 CONTINUE 
C 

IF (MDFLTL .EG* 0) GO TO 330 

EO 38A KFLOT «= JMX 

IF <MDFLOT< KFLOT) .EG. 0> GO TO 32A 

dPLOT = 0 

DO 308 J123 = Ij ii 

IF tNPLOT -EG. 1) MT1TLCJ123) = KTnLKJ123> 

IF <NFLOT .EQ. 2) MTITL<dI23> = MTI TL2< J 1 S3> 

IF <NFLOT *EQ* 3> MTITLCJ123) = KTITL3(J123> 

328 CONTINUE 

C 

DO 326 d = !* KUM 

IF <<UFLOTtNFLOT^d) *LT* YLOMD) *0F* (UFLOTCNFLOT* d) 

1 ,.GT* YHIMD)) GO TO 326 

JFLOT = JPLOT + 1 
DUKKYI CdFLGT) = TFLCKd) 

DU<KYY(dPLOT) = UFLOTCNFLCT^ J> 

326 CONTINUE 

IF CdPLOT *EG. 0) GO TO 32A 
PLOT AMPLITUDES OF DIFFEREN'T MODES. 

CALL 6RAFHS<IEUF/ 3000^ A^dPLOT# 51^ I TI CMD/ TMAX^YHIKD^ TKIK^ 

1 YLOKD^ ITT^MTITL, 21^20,DUKMYT* DUMKYY* 2*0* YLAFKL^ TITLE) 

324 CONTINUE 

IF <MDFL0T<4) *EQ. 0) GO TO 330 
JPLOT *= 0 
DO 332 J if NUM 

IF <<PFilT<<J) »LT* YLOMD) *0R« <PRITCJ) *GT* YHIMD)) GO TO 330 
JPLOT *= JPLOT + 1 
DWMYTC JPLOT) = TFLOTCJ) 

DWKYYCJPLOT) = PRiT'(J) 

332 CONTINUE 

IF (JFLOT .EQ. 0) 60 TO 330 
PLOT PRESSURE AMPLITUDE OF IT MODE. 

CALL GRAPH 5 ( IBUF*3000* 4/ JPLOT;> 51^ I TI CML» TMAX*YHIML* TWIN* 

1 YLOMD> I TT^ PRTl TL 2 29 DUMMY T:. DUKMYY ^2.0^ YLAEM D> TT TL E) 

330 CONTINUE 

REASSIGN PLOTTING PARAMETERS FOR NEXT SET OF PLOTS. 

JFLO T »= 0 
TMIN ^ TMAX 
TMAX *= TMAX + TLEL 

**«*++:ic*****;^ time HISTORY PRINTED OUTPUT SECTION **-:^=(c*+* *#+>!=**«♦* 
156 WRITE <6^6011) NSTEP^ TIME# (PRESStJ)# J >= 1#6)# 
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1 tAXVa-CJJ* J c 1,3)^ YPHI 

UINE s LINE + 1 

157 IF CTINE .GT. TQUIT) GO TO 850 
IF (LINE *LT. 52) 60 TO 155 

WRITE (6^6013) 

WRITE (6>6088) (ANGLE(J>, d = 1*6># CANGLE(d)» d *= 1*3) 

LINE c n 
C 

155 I e I + I 

IF (I .LT« LAST) GO TO 105 
C 

C *****:f:+*5|c+*)|t* LIMIT-CYCLE SECTION 

C 

C TEST FOR LIMIT CYCLE. 

K «s K + 1 

IF CCNTEST .EQ. 0) *0R. (MAXNO .LT. 80)) GO. TO 190 

UTOT =0.0 

DO 180 d = 0* 3 

JMAX •» MAXNO - d 

UTOT «= UTOT + ABS<UMAXCJKAX) ) 

180 CONTINUE 

UAUGCK) = UTOT/A.O 
IF (K .EQ. 1) GO TO 190 
CHANGE = UAUGCK) - UAUGCK- 1) 

ABSCHG = ABSCCHAN6E/UAUGCK) ) 

IF C ABSCHG .GT. ERR) GO TO 190 

TM = TIME/8.0 

ITM = TM 

ITM = 8+ITM ♦ 2 

TM = ITM 

TSTART = TM + TSTART 
TQUIT = TK + TQUIT 
TMIN = TSTART 
TMAX = TSTART + TDEL 
NTEST = 0 
C 

C RE-ASSIGN ARRAYS. 

190 DO 200 I = 1* NFl 

ILAST = LAST - NPl + I 
PRSCI) = PRSC ILAST) 

TICI) = TIC ILAST) 

DO 800 J = 1* NU 
UCI*d) = UCILAST*d) 

800 CONTINUE 
GO TO 100 
C 
C 

C *******♦+*>*:♦+ PRESSURE MAXIMA AND MINIMA PRINTOUT ♦*♦»**♦**#****♦* 

■C 

-250 WRITE <6*6083) ECNLOC)* ANGLECKLOO* MAXP 
LINE = 4 

DO 255 dST *= 1* MAXF* 8 
■ dSTABT = dST 
dSTOP = dST + 7 

IF CdSTOP »CT. MAXP) dSTOP = MAXP 
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S55 


V/RITt <6>60S4) (PKAX<J5# 0 = JSTART* OSTOFJ 

WRITR (6^6024) (TIMAXCJ)* 0 = dSTARl^ JSTOP) 

WRITJ (6^60143 

LINL LINE + 3 

IF CLINE »LT* 52) GO TO 255 

LINE = 0 

Vr'RITL C6>6013) 

CONTINUE ^ 

IF CCKOUT "EQ* 0) •0R» CNOKIT •££!• l>) GO TO 5 


*♦*♦ + *♦* 


PRESSURE MAXIMA PLOTTING SECTION *^*^*************** 


LE.TEPMINF. LARGEST VALUE OF PMAX-. 
AMFMAK s= 0«0 
DO 260 d = 1/ MAXP 

IF CFEAXCd) .LT» AMFMAX) GO TO 260 
AKFEAX s= FMAXCJ) 

£60 CONTINUE 


RANGE OF PLOT AND COORDINATE LABELING* 
ITK = AMFMAX + I-O 
AMPMAX »= ITK 

ITK = 1*0 + TIMAXCMAXP) /50*0 
TMAX - ITM + 50 
DELX = TMAX/10-0 
DELT = AMFMAX/ 10*0 


C 

C 


c 

c 


eliminate NEGATIVE VALUES* 
JFLOT 0 
DO 262 U = 1' MAXF 
IF <FMAXCJ>3 262j 264* 264 
264 OPLOT = JPLCT + 1 

DUMMYTC JPL0T3 = TIMAXCJ3 
DUMMYYC JPLCT3 = FMAXCJ3 
262 CONTINUE 


C^L GRAPHSCIBUF* 3000* 4*UpL0T* 101* 101* TKAX* AMPMAX* 0*0* 0*0* 
ITT* I TP* 2 1 * 14* I'LMjMYT* DUMMYY * DELX* ILLY » T1 TLE) 


1 


C 


GO TO 5 


C TUPN off PLOTTING ROUHNE* 

300 IF CNPT *EQ* CALL SHPAPG 


C 

c 

c 


^:^i^^Hc*****‘*** BEAD FOFMAT SPECI 


5000 FOFiMAT 

5001 FORMAT 

5002 FORMAT 

5003 FORMAT 

5004 FORMAT 

5005 FOFMAT 

5006 FORMAT 

5007 FOI<MAT 


C 12A63 

< ziFlO*0*3I5> 

< 51 5* 2F10.5* IX* A4) 
C 21 5) 

C2I 5*F15*6) 

<31 5*F15-63 

< 5F10.03 
(15* 2F10*0) 


FICATI ONS 


:4: *+:):*#* s|t ^ **#♦♦**** ^ ^ 
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5008 FOPKAT <71 S) ’ 

5009 FORMAT <7F10.0) 

5010 FORMAT <I5jAF10.5) 

5011 FORMAT <1S>2F10*5> 
501S FORMAT <F10«0> 

5014 FORMAT (AI5) 

5015 FORMAT <SF10*0*I5) 
5099 FORMAT <SF10»0) 


C 

C 

C 


«**xc4<***Ke«**:^ VRITE FORMAT SPEC! FI CATI OKS hp^****.*****************- 


6001 FORMAT 
1 

2 

6002 FORMAT 

1 

6003 FORMAT 
F0FF3AT 
FORMAT 
FORMAT 
FORMAT 


‘ ^F5 
.F5.2. 


I 


N NS# 7X# 3HSMN# 3X# 


6004 

6005 

6006 

6007 

6008 


C<#I 1# lOH# NJ#NF> 
-FIO. 5> 

DCNJ>NP#NQ) IS#I5/) 
#Fl0-5> 

INTERACTION INDEX = 


IS# I 5/) 


»F7.5# 


6009 

6010 
6011 
6012 

6013 

6014 

6015 

6016 


<IH1#9H GAMMA = # F5« 3# 5X# SHOE 
5X#5HZE = , FS. 5# 5X#8HZC0ME *= 

5X#8HNJMAX' = #I2//> 

<2X#29HNAME J L M 
7HJM< $MN>/> 

<2X# A4# 51 5# 2F10-5) 

(1H0#26H NOMEER OF COEFFICIENTS 
< 8X# 2HC< # 1 1# IH## 1 2# IH# #T2# 4H5 = 

( IHO# 38H NUMBER OF COEFFICIENTS 
C2X#2HD(#I2#1H##I8#1H##I2#4H> < 

FORMAT< IHI#45H COMBUSTION FARAMETERSt 

12X# 1IHTIME-LA6 e # F7 . 5/2X# IThMOTOR PARAMETERS: # 19X# 

I BHGAMMA *= #F7«5#E3H EXIT MACH NUMBER = »F7.5# 

I 22H LENGTH/DIAMETER #F7.5//5 

C2X# 18HINITIAL CONDITIONS//) 

C IHO# 5X# 1HJ#8X# EHYR#8X#8HYI#7X#3HEPS# 7X# 3HETA//3 
C2X#I5#F18*5# lOFlO-5) 

1 1H0> 

ClHl) 

CIH ) 

C2X#I5#4F10*5) 

<1H1#36H INITIAL CONDITIONS ARE OF THE FORK:// 
8X#49HU<1#J) AC< J)*COS<FREG*T) + ASCO)*S1N<FREO* T) )# 
14H * EXP<DAMP*T)///6X# 1HJ#8X#7HDAMFIK6# 

6X#9HFREQUENCY# lOX# 5HACC J># lOX# SHASCJ)//) 

<8X#I5#4F15«e/) ,, 

<1HI#46H COEFFICIENTS FOR CONFUTATION OF VALL PRESSURE# 
lOH VAVEF0EM5///43X#27HC0EFFICIENT$ IN SERIES FOR:// 

82X# 5HTHETA# 10X#Z|HTIME# 1 OX# 5H THETA# 10X#5H AXIAL/ 

6X# 1HJ#9X#-IHZ# 3X#9H< DEGREES)# 5X# lOHDERIVATI VE# 
5X# lOHDERlUATIVE# 5X# lOHDEhlUATl VE//) 

(8X#I 5#F10»3#F18. 1#3F 15.7> 

(26X# 17HINJECTQR PRESSURE# 1 4X# 1 5HN0ZZLE PRESSURE# 

I2X# 21HN0ZZLE AXIAL VELOCI lY/SX# 4HSTEF# 8X# 4HT1ME# 
F5*0#5H deg# # F5* 0# 5H DEG»»F5»0#5H DEG## 

F5#0#5h DEG##F5#0#5H DEG##F5#0#5H DEG## 

F5#0#5H DEG.#F5*0#5H DEG##F5#0#5H DEG## 6X# 4HYFHI //) 
<1H1#38H PRESSURE MAXIMA AND MINIMA AT: t #F5*2# 

IIH THETA #F4*1/19H VALUES COMPUTED* #I3//> 

(IH #7X#8F13*6) „ ^ . 

C2X//8X# 37HTHE TRANSIENT BEHAVIOR IS CALCULATED*) 

< 2X//2X# 39HTHE LIMIT- CYCLE BEHAVIOR IS CALCULATED*) 
(8X//2X# 33HTHIS RUN PRODUCES PLOTTED OUTPUT*) 
C8X//2X#"THE PHANTOM ZONES ARE ELIMINATED*") 


FOFiMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 


6017 FORMAT 

6020 FORMAT 
1 

2 

3 

4 

6021 FORMAT 

6022 FORMAT 
1 

2 

3 

4 

6023 FORMAT 

1 

6024 FORMAT 

6025 FORMAT 

6026 FORMAT 

6027 FORMAT 

6028 FORMAT 
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6030 FOFnKAT ( 2X> ‘T-EOFLET MOMENTUM SOUECE IS NEGLECTED”/) 

6031 FOFINAT ( SX^ ”DEOFL ET MOMENTUM SOURCE IS INCLUDED”/) 

6032 FORMAT <2X* "NOZZLE NONLINEAPI TI ES NEGLECTED"/) 

6033 FORNAX <2X, "NOZZLE NONLINEARITIES INCLUDED"/) 

603A FORMAT < IHO# 8X> I H J, 1 OX^ 2HGR> 1 OX, SHGI // ) 

6035 EOSMAT ( 5X, I 5, 2F 1 2. 5) 

END 
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SUBROUTINE PHICFSCNFjZ* THETA^CTj CTH* CZ ) 


This SUBROUTINE COMPUTES THE COEFFICIENTS NEEDED TO 
CALCULATE THE WALL PRESSURE PERTURBATION. 

NP IS THE INDEX OF THE COMPLEX SERIES TEI30. 

Z IS THE AXIAL LOCATI ON<i 

THETA IS THE AZIMUTHAL LOCATION. 

CT IS THE COEFFICIENT IN THE SERIES FOR THE TIME DERI UATI UE 0 
THE VELOCITY POTENTIAL. 

GTH IS THE coefficient IN THE SERIES FOR THE THETA DERIVATIVE 
OF THE VELOCITY POTENTIAL. 

CZ IS THE COEFFICIENT IN THE SERIES FOR THE AXIAL DERIVATIVE 
OF THE VELOCITY POTENTIAL. 


C 


C 


c 


PARAMETER 

COMPLEX 

1 

COMMON 


MX = 5 

Cl. CZ. CAXI. CAXIZ. CRAD. CAZI. CAZI TH. 
BtKX). CT^ CTH> CZ 

/BLK2/ MCMX). NSCMX>^ SJCMX>. B 


Cl = (0.0. 1.0> 

CZ = CMPLXCZ.O.O 

CAXr «= CCOSHCCl ♦ BCNF> * CZ5 

CAXIZ = Cl * BCNF) ♦ CSINHCCI * bCNP) ♦ CZ> 

CHAD = CMFLX(SJCNF>.0-0> 

EM e MCNF) 

ARB J= EM + THETA 
FSIN = SINCARG) 

FCOS *= C0S(ARG5 
AZI p FCOS 

IF (NSCNP> .EG. 1> AZI = FSIN 
AZITH c EM ♦ FCOS 

IF CNSCNP) .EQ. 2) AZITH » -EM * FSIN 
CAZI «= CMPLXCAZI.0.05 
CAZITH c CKPLXCAZITH. O.OJ 

CT = CAZI ♦ CAXI * CRAD 
CTH = CAZITH * CAXI * CHAD 
CZ s= CAZI ♦ CAXIZ * CRAD 

RETURN 

END 
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SUBROUTINE PESVEL.C UEAR> VTH> V2 3 


C 

C THIS SUEFOUTINE COMPUTES THE WALL PHESSLKE AND VELOCITY. 

C 

C UEAE IS THE LOCAL AXIAL STEADY STATE MACH NtiMBER. 

C UMS IS THE DEEIUATIUE OF THE MACH NUMBER EOR THE CASE 

C WHEE! DROPLET M0MENTLa>1 SOURCES ARE INCLUDED. 

C Y IS THE ARRAY CONTAINING VALUES OF THE MODE- AMPLI TUDE 

C FUNCTIONS AND THEIR DERIVATIVES. 

C PIS THE VALUE OF THE WALL PRESSURE PERTURBATION. 

C VTH IS THE TANGENTIAL COKPONEN.T OF VELOCITY AT THE WALL. 

C VZ IS THE AXIAL COMPONENT OF VELOCITY AT THE WALL. 

C 

PARAMETER MX2=10# MXA=20 

DIMENSION YCKXAJ* SUMCAJ> SUKS0C3) 

COMMON /-BLK3/ NJKAXjp NLMAX> GAMMA. COEFC3.KXSJ 

C 

DO 10 I 1. A 
SUM<I> =0.0 
10 CONTINUE 
C 

bo SO I = 1. A 
DO 20 J = 1 j KJMAX 
JY = J 

IF a .£Q. i) JY = J + NJMAX 
II = I 

IF (I .EQ. A) II = 1 

SUMCI) = SUMCI) + Y<dY> ♦ COEF(lI.d) 

20 CONTINUE 
C 

FLlK = SUM<1) + UEAR*SUMC35 + UMS=*=5UK<A> 

PNL = 0.0 

IF CNLMjAX -EO. 0> go to AO 

DO 30 I = 1. 3 

SUKSOCI) = SUM(I) St SUM(I> 

30 CONTINUE 

PNL = 0.5 + (£UM£0<2> + SUMSG<3) - SUMSOC 1) > 

C 

AO P = -GAM.MA + (PLIN + PNL) 

VTH = SUM<2) 

VZ « SUMC3) 

C 

RETURN 

END 
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SUBROUTINE RHSCNU# 1 1# U* UP5 

PARAMETER MX«=5^ MX2=10j MX4«S0> «X2£G=100 

DIMENSION UCNU5». UFCNU5 

COMMON RV<MXSj. 4>» ’.CC/i^KXa^MX’e’j^ D<MX2>MX2Se)* 

1 KFMAXC4#MX25* I C< 4>MXe-MX2)# KPQMAXtKXPI > 

2 . IDPCMX2»MX2SQ)* 1 LQCMXB»MX2£G3 

COMMON /ELK3/ NJMAX* NLMAX* GAMMA, C0EF<3';MX2 

COMMON /NLTEFM/ NOZNLS> £XTRACMX2,A> 

DO 10 NJ = 1, NJMAX 
KJP c NJ + NJMAX 
UF<NJ> = UCNJP) 

sLi = o;o 

. SL2 =0*0 
SL3'= 0.0- 
£LA = 0.0 
SNL = 0.0 
MAX c KFMAXCI>NJ> 

IF (MAX .EQ. 0) GO TO 8S 
DO 20 KP = 1, MAX 
NF c lCCi>NJjKP) 

SLI «= SLI + <CC1,NJ,KF5 * U<Np)J 
20 CONTINUE 
25 MAX s KR1AXC2 j»NJ3 

IF CMAX -EO. 05 GO TO 35 

DO 30 K-P «r 1, KAX 

NPP c ICCe#NJ*KF5 '+ NJMAX 

SL2 = SL2 + CCC2,NJ,KF5 * U<NFF5> 

30 CONTINUE 

35 MAX KFMAXC3»NJ5 

IF <MAX‘-.EQ. 05 GO TO A5 
' ’ DO AO-KP *= 1> MAX 

NF B IC<3>NJ,KP5 

ELS = SL3 + fCC3»NJ,KP5 * RV<NP,I15 5 
AO CONTINUE- . - 

A5 IF CN0ZNL2 .NE. 15 GO TO 65 
MAX = KFMAX(A»NJ5 
'■ IF CMAX .EQ. 05 GO TO 65 
DO 60 KP B MAX - 
NP = ICCA;iN0/KF5 

SLA = SLA + CCCA,NJ,KP5 ♦ EXTRACNF, 1 1 5 5 
60 CONTINUE 

65 IF CNLMAX -EG. 03 GO TO 55 
MAX B KFQMAXCNJ> 

IF CMAX -EG. 03 GO TO 55 
DO 50 KFQ B if MAX 

NP = IDPCNJ,KFO) . . 

NQP = IDOCNJ,KPG> + NJMAX 

SNL = SNL + CDCNJ,KFQ> * U(KF> * U<NeP>3 
50 CONTINUE 

55 UPCNJP5 = “CSLl + SL2 + SL3- + SLA + SNL3 
■ 10 CONTINUE 
RETURN 
END 
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COMPILER <FLD=ABS) 

SUBRO UTI N E G RAPHS C I BUF> KLOC^ LDEV ^ N TOT^ NTI CK ^ CY ^ 

1 X W AX':» Y K AX , XM I K YM I K , I T 1 TL X ^ I T1 TLY , L TI TL / L TI TLY ^ X A HR AY > 

2 YABBAY>DELXjDELY>TITLE'> 


c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IDENTIFIER 


meaning 


IBUFJ 

NLDC: 

LDEU: 

NTOTs 

NTICX: 

NTICYJ 

XMAX! 


ADDRESS OF BUFFER AREA FOB PLOT OUTPUT 
NUKEEB OF LOCATIONS IN BUFFER AREA (>=2000> 
LOGICAL DEVICE NUMBER FOR PLOT 
NUMBER OF POINTS TO BE PLOTTED 
NUMBER OF TIC MARKS ON ABSCISSA <>=2) 
NUMBER OF TIC MARKS ON ORDINATE 0=2) 

UPPER LIMIT OF ABSCISSA DOMAIN 


C YMAX: 

C XMIN; 

C YMINi 
C ITITLX: 
C ITITLY: 


UPPER LIMIT OF ORDINATE RANGE 
LONER LIMIT OF ABSCISSA DOMAIN 
LOWER LIMIT OF ORDINATE RANGE 
ABSCISSA LABEL 
ORDINATE LABEL 


1YPE 


INTEGER 

INTEGER 

INTEGER 

INTEGER' 

INTEGER 

INTEGER 

REAL 

REAL 

FsEAL 

REAL 

FIELDATA ARRAY 
FI EL DATA ARRAY 


C LTITLX 
C LTITLY 
C XARRAY 
C y ARRAY 
C DELX: 

C 

C DELYi 
C 

C TITLE: 
C 

C 


NUMBER OF CHARACTERS IN ITIILX 
NUMBER OF CHARACTERS IN ITIILY 

ABSCISSA POINTS IK TERMS OF XMIN-XMAX COORD'S 
ORDINATE POINTS IN TERMS OF YMIN-YMAX COORD'S 
INTERVALS OF ABSCISSA TIC MARK LABELING 
IN TERMS OF XMIN-XMAX COORDINATES 
INTERVALS OF OEBINAIE TIC MARK LABELING 
IN TERMS OF YKIN-YMAX COORDINATES 
LABEL FOR THE WHOLE. RUN 


INTEGER 
INTEGER 
REAL ARRAY 
REAL ARRAY 

REAL 

REAL 

FI HLDATA ARRAY 


DIMENSION IBUFCKLCC) jXARRAVCNTOT)*YARBAYCNTOT).. ni ILXC I )# 
1 I TITLY C DjYDITC 100) 

DIMENSION TITLECl) 


FIXED BASIC PARAMETERS 


LOGICAL ZERO 

DEFINEZERC'=NDEC*LT*G.AND-ABS<FPN) .LT . *5 
1 *0R*NDEC.GT*0»AND.AB5CFPN) .LT<>5.=^' 10 *+=*=C -NLEC- 1) 
DEFINE DNDEC=Ni)EC-FLD(0.» 36^ ZEfiO)*NDEC-FLD< 0.» 36>ZER0) 
DEFINE IFIX<FABG) = 1NT<FAP.G+- 5) 

DATA J/1/ 

HEIGHT/. 105/ 

IN TEC/ 1/ 


DATA 

DATA 

DATA 

DATA 

DATA 


AHSCIS/8./ 
OFtDINA/6./ 
I CODE/- 1/ 






QU, 





DATA TOPMAR/1,/ 
data B0TMAR/1*5/ 

REAL LEFMAR 
DATA LEFM AR/ 1.9/ 

DATA RYTMAR/1.1/ 

DATA FACT/1./ 

DATA MAXIS/1/ 

DATA MLINE/1/ 

DATA HTLAB/.105/ 



c 

c 19 INITIAL COMPUTATION OF DERIVED PARAMETERS 

C AND INITIAL PLOTS CALL 

0 20 SKIPS PRELIMINARIES FOR 2ND AND SUBSEQUENT CALLS 
C 



GO TO <19^20)>tJ 
19 YDIT(1> « 3./19. 

TICKLE a. HEIGHT/2. 

EOTFAC = - 3*/lA. + HEIGHT - 4./7« * HEIGHT 

STARTL e 6 * height + ROTFAC + TICKLE 

SEPLAB s: STARTL + 1*5 * HEIGHT 

SYMBLH = 0.070 

REAL LAB SEP 

LABSEP = 4. ^ HEIGHT 

ASTART = 2. * HEIGHT 

DO 1 I 2^ 100 

1 YDITCn = YDITCI - 1) + C2 + M0D<1,21 + 15 / 19 . 

VDITC1005 = YDITClOO) + .5 

CALL PLOTSCIBUF^NLOC>LDEV> 

CALL FACTORCl.) 

0=2 

SYMBOL <HE1GHT^36 * HEIGHT + 5. 5^KEIGHT> TITLE^ 270.>72> 
CALL# PLOTCl^^ “ *5^ 35 

3 DO 2 I K I, 100 

e CALL PL0T<0. ^YDITCI 3 - K0DC1#235 

DO 33 I = 1> 100 

33 YDIT<1) = YDITCI) - ABSCI S - RYTMAR 


c RESET ORIGIN 

C 



XPAGE = BOTMAR + ORDINA 
GO TO 2019 

20 XPAGE = BOTMAR + ORDINA -i- TOPMAR 
2019 CALL VHERE<RXPAGE^RYPAGE.FACT> 

YPAGE = HYPAGE - LEF1*;aH - 
CALL PLOT<XPAG*Ej.YPAGE, - 3> 

CALL FACTOR (FACT) 
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Q 

c 

C DRAW AXES AND LABELING MAXIS TIMES 

C 

C 

DO 100 I = l^MAXIS 
100 CALL MY AXIS 

C 

C 

C DRAW POINTS^ OPTIONAL CENTERLINE^ AND PAGE SCI SSOEL 

C MLINE TIMES 

C - . 

c 

BO £0O 1 = 1>MLINE 
eoo CALL MYLINE 

RETURN 

C 

C 

C ENTRY POINT SHPARG 

C TERMINATE PLOTTING SEQUENCE 

C 

C 

ENTRY SHPARG 

CALL UHERE<KXPAGE^KYPAGE^ 1) 

CALL PL0TCKXPAGE>RYPAGE^999> 

RETURN 

C 

C SUBROUTINE MYAXIS UNTERNALJ 

C 

c 

SUBROUTINE MYAXIS 

STARTL = 6 HEIGHT + ROTFAC + TICKLE 
IMAX = IFIXICYKAX - YMIN)/D£LY5 
TICSEP = ORDINA/<ABS<NTICY) - 1) 

CALL DENr>EC<YMAX;.DELY#NDEC) 

K = 1 

N = < ABS<NTICY)/IMAX) - 1 + MODI ABSCNTI CY ) ^ 2) 

DO 9 I = 0,IMAX 
GO TO Cn^l2>*K 

11 IFC2 * I.LT.IMAX5G0 TO 12 

CALL AXLABCO* ^ I TI TLY ^ LTI TLYj HTLAB) 

K = £ 

12 FPN = YMAX - I + DELY 

IFCZEROJFFN = 0- 
TMID = 1. 

XPA6E = - I * ORDINA/IMAX - .5 * HEIGHT 

IF( FPN) 1 13^ 122^ 1 18 

113 IFCNDEC - 2)115^114^112 

114 YPAGE = STARTL@5CHAR 
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60 TO lia 

US IFCNDEC - nU7*U6*U2 

116 YPAGE = STARTL ~ HEIGHTS ACHAR 

GO TO 112 

117 IFCABSCFPNJ - 100. > 119» 1 16^ 1 16 

119 IFCABS<FPN> - 10 - ) 1 20^ 121^ 1 21 

120 YPAGE = STARTL - 3 * HEIGHTS2CHAR 

GO TO 112 ■ !- 

121 YPAGE = STARTL - 2 * HEIGHT@3CHAR 

GO TO 112 - ’ . j 

122 YPAGE = STARTL - A * HEIGHTS ICHAR 
GO TO 112 

118 IFtNDEC - 21123^116.112 

123 IFtNDEC “ 1 J125> 12A. 112 

124 IFCFFM - 10.1121. 116. 116 

125 IFCFPN - 10-1122.120.126 

126 IFCFPN - 100-1120.121.127 

127 IFCFPN - 1000-1121.116.128 

125 IFCFPN- 10000.1116.114.114 

112 NNDEC - DNDEC ’ 

CALL NUMBER CXPAGE.YPAGE.HEIGHT.FPN.270-.KNDEC) 

XPAGE = - I * CORDINAVIMAXl- 

DO 10 JJ - l.N 

YPAGE = TICKLE * TMID 

CALL PLOTCXPAGE. YPAGE. 31 

YPAGE = YPAGE + C - 1 + 1/IWAX * .51 

CALL PLOTCXPAGE. YPAGE. 21 

I FC I /I MAXI 110. 1 10.9 

110 YPAGE = 0 

CALL PLOTC-XPAGE.YPAGE. 31 
XPAGE = XPAGE - TICSEP 
CALL PLOTCXPAGE. YPAGE. 21 
TMID = .5 

10 CONTINUE 

9 CONTINUE 

K » 1 

IMAX s= IFIXCCXMAX - XMINl/DELXl 
TICSEP = ABSCIS/CKTICX - 11 
XPAGE = >- ASTART - ORDINA 

CALL DEMDECCXMAX.DELX.NDEC) 

' ' DO 25 I = O.IMAX' ■' ■■ 

• STARTL = - I * ABSCIS/IMAX 

GO TO C24. 251.K 

24 IFC2 * I-LT.IMAXIGO TO 25 

■ CALL AXLABC270.;ITITLX.LT1TLX.HTLAB> 

K = 2 

XPAGE = “ ASTART - ORDINA 

25 FPN '= XMIN + I + DELX 
IFCZEROIFPM = 0. 

IFCFFN)813.822^818 ' •' ' 

813 IFCNDEC - 218T'5.814.23 

814 YPAGE ■ STAHTL '-i-'- 1'6'.V7. - * 

GO TO ,23 

815 ••'I'FCNDEC - 


1)617.616.83 



816 

817 

819 

820 

891 

822 

816 

823 

824 

825 

826 

827 

828 
23 
28 


111 


27 

26 

C 

C 

C 

C 

C - 


17 


YPA6E = STARTL + 25*/ 14. * HE1GHT©4CHAR 
GO TO 23 

IF(ABSCFPN) “ 100. 1819, B16>816 
IECABS<FPN) - 100820^821^821'' 

YPAGE = STARTL + 11./ 14. =(= ' HEIGHTe2CHAR 
GO TO 83 

YPAGE = STARTL + 9./7- * HEIGHTOaCHAR 
GO TO 23 

YPAGE = STARTL + 2./7* ♦ HEIGHT© ICHAR 
GO TO 23 

IFCNDEC - 2>823>816j23 
IFCKIDEC - 1)825^824^23 
IF<FPN - 10.)82U816^816 
IFCFFN - 10. >822^820^826 
IFCFPN - 100.>820^821>627 
IFCFPN - 1000. >821j 816^828 
IFCFFN - 10000. >816j814j814 
NNDEC = DNDEC, 

CALL NOKBERCXPAGE^YPAGE^HEIGHT^FPN^ 270.>NNDEO 
N = CNTICX/IMAX3 - 1 + l^ODCNTI CX^ 2> 

BO 26 I = 1KAX>0j - 1 
TMID = 1. 

YPAGE = - 1 + ABSCIS/IKAX. 

BO 27 OJ = l^N 

XPAGE = - ORDINA - TICKLE + TMID 

CALL. PL0TCXPAGE>YPAGE>3> 

XPAGE = XPAGE + < TICKLE + FLDCO^ 36^ I .NE*0> * TICKLED * TMID 
CALL PL0T(>:PA6E>YPA6E>2> 

IFC1>111^26^ 111 

XPAGE = - ORDINA 

CALL PL0TCXPAGE;.YPAGE^3> 

YPAGE = YPAGE + TICSEP 

CALL PLOTCXPAGEjYPAGE# 23 

TMID = .5 

CONTINUE 

CONTINUE 

RETURN 


SUBROUTINE MYLINE CINTEHNAL) 


.53/ U » * 99.3 


SUBROUTINE MYLINE 
I TOP = IFIXCCABSCIS ^ RYTMAR + 

IBOT = I Fixe RYTMAR/ 5 1 . 99 . 3 

DO 17 I = 1/NTOT 

XPAGE == (YARRAYC13 - YMAX3/CYMAX - YMIM) 
YPAGE = CXMIM - XARHAY C 1 3 3/ (XMAX - lU^IW) 
CALL SYMBOLCXPAGE^YPAGE:. SYMELH# INTEQ> 270 
IFCNT1CV.GE.03G0 TO 22 
XPAGE = - ORDINA/ 2. 

YI>AG£ = - ABSCIS 

CALL PLOTCXPAGE^YPAGE^ 33 
DO 18 I = IBOT^ITOP 


ORDINA 
* ABSCI S 
.^ICODE) 


of B 

■ QUAUry 
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18 CALL PLOTCXPAGE^YDITCIl^S - M0DCI>2)) 

22 XPAGE = TOFMAR 

YPAGE = - ABSCIS - EYTMAR - .S 

CALL PLOTCXPAGE^YPAGEj 3> 

DO 21 I = 1^100 

21 CALL PL0T<XPAGE>YDITCI1,3 - MODCl'^2)? 

RETURN 

C 

C 

C SUBROUTINE AXLAB (INTERNAL) - 

C 



SUBROUTINE AXLABC ANGLE^ I BCD^ MCHARX>HEI GHT) 

DIMENSION 1BCDC7) ’ 

LOGICAL S' 

INTEGER QSQ/' S'/ 

K = 2 

NCHAR = NCHARX 
S = *FAL£E. 

IF(ABSCAK6LE) .GT.. DGO TO 30 
XPAGE = - ORDlNA/2. ~ NCHAR * HEIGHT/2 

yPAGE = SEPLAB 
GO TO al- 
so XPAGE' = - OBDINA - LABSEP 

YPAGE = - ABSCIS/2. + NCHAR HElGHT/2 

31 LSTART = 6 sK M0D(NCHAR^6) ~ 12 
IF(LSTABT.EQ. - 12)LSTART = 2A 
LOOK = NCHAR/6 + 1.1 
IFCLSTART.EG. *• 6) GO TO 13 

1F(FLD<0^ 12;* S' ).EQ.FLDCLSTART* 12^ IBCDCLOOK)))GO TO 15 
GO TO 14 

13 IFCFLD(0^6j ■ j • ) .NE.FLD(30^6> IBCDCLOOK - i)))GO TO 14 
JF<FLDCO^ 6^ • S' ) .NE.FLDCO# 6 j IBCDCLOOK))>GO TO 14 

15 NCHAR = NCHAR - 1 
S = .TRUE-. 

14 CALL SyMB0L(XPA6E#YPAGE*KEIGHT^ IBCD>ANGLE>NCHAR) 
IFCSICALL SYMB0LC999.^999.^ 2 * HEI GHT/3> QSQ* ANGLE, 2) 
RETURN 

C 

C 

C SUBROUTINE DENDEC (INTERNAL) 

C 

C 

SUBROUTINE DENDEC ( GMAX* DELO>WDEC) 

IF(INT( ABS(QWAX) ) .GE* 10)G0 TO 5 
IF(AMOD(ABS(QMAX - DELQ) * . 1 ) . GE. . 0 1 )G0' TO 7 
NDEC = 1 
r/jrtSjRETURN 
5 NDEC = - 1 

RETURN 

7 NDEC = 2 

RETURN 

END • - - 
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